Моделирование методом Монте-Карло

Содержание

  • Откуда берутся доходы?
  • Хорошо ли это смотрится на мне?
  • Куда вы идете с этим?
  • Дойная корова Шредингера
  • Оптимизация параметров и анализ чувствительности для оптимизации качества процесса в Python
  • Прогнозирование будущей стоимости собственного капитала с помощью логнормального распределения
  • Построение вероятностной оценки риска с использованием моделирования по методу Монте-Карло с помощью Python MCerp
  • Моделирование по методу Монте-Карло и треугольное распределение в Python
  • Моделирование по методу Монте-Карло и равномерное распределение в Python
  • Метод Монте-Карло: ретроспективное моделирование
  • Игровая площадка квантов
  • Имитационное моделирование методом Монте-Карло

Часть 1: Откуда берутся доходы?

Как распределяются возвраты? Как мы можем откалибровать наши ожидания в отношении будущих событий? Это сложный вопрос, на который многие пытались ответить, с противоречивыми ответами и с разрушительным эффектом. В этой серии статей мы попытаемся ответить на этот вопрос, используя оценку максимального правдоподобия, правило Байеса и другие статистические методы. Для статистического анализа мы будем использовать библиотеки python, в частности модуль Scipy.stats, API Stats Models, а также основные Numpy и Pandas.

Обо всем по порядку

Мы загрузим и изучим набор данных ценовой истории для тикера на Египетской бирже, в этом примере мы будем использовать ETEL:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

ETEL = pd.read_csv('myfolder/ETEL Historical Data.csv')

display(ETEL)

Сначала мы должны отформатировать наши данные, чтобы мы могли их проанализировать, для этого мы устанавливаем дату в качестве индекса DataFrame, нам также нужно изменить порядок данных так, чтобы самая ранняя дата была первым элементом:

def clean_stock(tickername):
tickername['Date']= pd.to_datetime(tickername['Date'], format='%m/%d/%Y')
tickername.sort_values('Date',inplace = True)
tickername = tickername.set_index('Date')
return tickername

ETEL_clean = clean_stock(ETEL)
plt.plot(ETEL_clean['Price'])

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

def log_returns(history):
return np.log(history/history.shift(1)).fillna(value = 0)

log_rets = log_returns(ETEL_clean['Price'])

Затем мы можем визуально рассмотреть его с помощью гистограммы:

plt.hist(log_rets, bins = 100, density = True)

Выглядит вполне нормально, правда? Примерно в форме колокола, ничего необычного, давайте проведем тест на нормальность на этих данных, а затем продолжим, чтобы оценить, насколько это близко к нормальности. Для нормальности используем тест Шапиро-Уилка:

from scipy.stats import shapiro

shapiro(log_rets)

Получаем следующий результат:

Тест возвращает p-значение, равное примерно 0, ниже 0,05 или любой значимый уровень значимости, поэтому у нас есть достаточно доказательств того, что эта выборка не происходит из нормального распределения.

Но почему так? Гистограмма ясно показывает, что это в какой-то степени нормально, верно?

Насколько это плохо?

Ключ лежит в выбросах слева и справа от гистограммы, чтобы понять, почему, мы строим график Q-Q данных, чтобы упростить понимание визуального элемента, мы стандартизируем данные:

import statsmodels.api as sm
# here we subtract the mean from the data and then divide by the standard deviation
# so that every return occurs at the corresponding standard deviation
log_rets_stddized = (log_rets - log_rets.mean())/log_rets.std()

sm.qqplot(log_rets_stddized, line ='45')

Если бы выборка действительно была из нормального распределения, мы ожидали бы более или менее идеальной корреляции между теоретическими и выборочными величинами. График Q-Q представляет собой сравнение того, сколько точек данных присутствует в каждом квантиле выборочных данных, с тем, сколько их должно быть, если выборка была совершенно нормальной. Область вокруг среднего, по-видимому, выстраивается довольно хорошо, однако по мере удаления наблюдается последовательный дрейф за порог двух стандартных отклонений. (Порог может быть произвольным, однако то же самое явление присутствует в большинстве акций по всему миру) Еще большую тревогу вызывает то, что доходность превышает 4, 6, почти 8 стандартных отклонений.

Это ПЛОХО, ПЛОХО

Реальная проблема здесь в том, что эти выбросы нельзя игнорировать ни в коем случае, требуется больше доказательств, чтобы опровергнуть нормальность, чем доказать ее, даже в самых широких приближениях, чтобы довести дело до конца, мы позаимствуем доказательство с помощью закона Байеса из книги Николаса Нассима Талеба «Статистические последствия жирных хвостов»

Мы оцениваем вероятность того, что событие произойдет на расстоянии от 2 до 5 стандартных отклонений, один раз по эмпирическим данным, и сравниваем ее с вероятностью из нормального распределения.

from scipy.stats import norm
def Bayes_Normality(Prior_probability, significance):
# we calculate the empirical probability of returns above a certain number of standard diviations
Probability_sigma = len(log_rets[log_rets > significance*log_rets.std()])/len(log_rets)
# then we calculate the theoretical probability for a normal distribution
Probability_sigma_normal_distribution = norm.sf(significance*log_rets.std(),loc = log_rets.mean(), scale = log_rets.std() )
# this is how sure we are of Normality
Prior_Prob_Gauss = Prior_probability
# this is how unsure we are of Normality
Prior_Prob_Not_Gauss = 1 - Prior_Prob_Gauss
# here we apply Bayes Rule in accordance with the book
Posterior_Prob_Gauss = Probability_sigma_normal_distribution*Prior_Prob_Gauss/(Probability_sigma_normal_distribution*Prior_Prob_Gauss + Probability_sigma*Prior_Prob_Not_Gauss)
return Posterior_Prob_Gauss

prior = 0.9
for sigma in range(2,6,1):
posterior = Bayes_Normality(prior, sigma)
print('Significance: {}, Posterior Probability: {}'.format(round(sigma,2), round(posterior,2)))

Даже если мы на 90% уверены в Нормальности, мы получаем следующие результаты:

Здесь мы видим, что, чем более вероятны экстремальные события, тем менее верно предположение о нормальности, настолько, что при 5 стандартных отклонениях практически невозможно, чтобы данные были нормальными. Обратите внимание, что, хотя события стандартного отклонения >5 крайне маловероятны при нормальности, они являются вполне реальной возможностью в реальной жизни, сравнение вероятностей показывает:

Это буквально в тысячу раз менее вероятно при предположении о нормальности.

Таким образом, мы доказали, что доходность не является нормальной, одна и та же серия тестов может быть выполнена практически с доходностью любой акции, большинство из которых дают один и тот же результат снова и снова. Теперь вопрос в том, если они ненормальны, то откуда берутся эти доходы? На этот вопрос нельзя однозначно ответить, однако мы можем предложить более вероятные альтернативы, такие как распределение Стьюдента, распределение Коши, а также стабильное распределение, которые будут рассмотрены в следующей статье.

Часть 2: Хорошо ли это смотрится на мне?

В предыдущей статье (здесь) мы привели несколько аргументов и объяснений ненормальности доходности акций. В этой статье мы стремимся найти альтернативное решение (решения) головоломки с использованием оценки максимального правдоподобия (MLE). Если вы не знакомы с MLE, не волнуйтесь! Все будет объяснено интуитивно и пошагово.

Максимальное правдоподобие что теперь?

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

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

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

Подгонка по вашему вкусу

Наша цель здесь состоит в том, чтобы найти наиболее подходящее распределение для данных, нашими кандидатами здесь являются распределение Стьюдента-t, стабильное распределение и распределение Коши. Для сравнения мы также попытаемся использовать нормальное распределение с единственной целью показать, насколько оно некомпетентно.

Первое, что нам понадобится, это функция, которая вычисляет логарифмическое правдоподобие с учетом данных и параметров, а именно:

К счастью для нас, нам не нужно вручную вычислять это для каждого дистрибутива, Scipy.stats приходит нам на помощь со встроенными функциями. Во-первых, мы импортируем наши данные:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

ETEL = pd.read_csv('myfolder/ETEL Historical Data.csv')

def log_returns(history):
return np.log(history/history.shift(1)).fillna(value = 0)

def clean_stock(tickername):
tickername['Date']= pd.to_datetime(tickername['Date'], format='%m/%d/%Y')
tickername.sort_values('Date',inplace = True)
tickername = tickername.set_index('Date')
return tickername

ETEL_clean = clean_stock(ETEL)
log_rets = log_returns(ETEL_clean['Price'])

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

from scipy.stats import norm, t, levy_stable, cauchy

def Gaussianlog_log_likelihood(params):
mu = params[0]
sigma = params[1]
data = log_rets
n = len(data)
dt = 1
log_likelihood = norm.logpdf(np.array(data),loc = mu, scale = sigma).sum()
return -log_likelihood

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

Для этой задачи минимизации мы будем использовать генетический алгоритм для решения, используя пакет python geneticalgorithm (получите его здесь, читайте об этом здесь и здесь):

from geneticalgorithm import geneticalgorithm as ga

varbound=np.array([[-0.1,0.1],[0,0.1]])

algorithm_param = {'max_num_iteration': 1000,\
'population_size':1000,\
'mutation_probability':0.2,\
'elit_ratio': 0.05,\
'crossover_probability': 0.5,\
'parents_portion': 0.3,\
'crossover_type':'uniform',\
'max_iteration_without_improv':None}

modelGaussian=ga(function=Gaussianlog_log_likelihood,
dimension=2,
variable_type='real',
variable_boundaries=varbound,
algorithm_parameters=algorithm_param)

modelGaussian.run()

Затем решатель дает нам:

Мы нашли наиболее подходящее среднее значение и стандартное отклонение с логарифмической вероятностью около 7800. Мы сохраним это в памяти для сравнения с другими дистрибутивами. Теперь мы повторяем описанный выше процесс с t-распределением Стьюдента:

def T_distribution(params):
  dof = params[0]
  loc = params[1]
  scale = params[2]
  data = log_rets
  Log_likelihood = np.sum(t.logpdf(np.array(data),df = dof, loc = loc, scale = scale))
  return -Log_likelihood

varbound_T= np.array([[0,10],[-1,1],[0,1]])

algorithm_param = {'max_num_iteration': 1000,\
                   'population_size':1000,\
                   'mutation_probability':0.2,\
                   'elit_ratio': 0.05,\
                   'crossover_probability': 0.5,\
                   'parents_portion': 0.3,\
                   'crossover_type':'uniform',\
                   'max_iteration_without_improv':None}

modelT=ga(function=T_distribution,
         dimension=3,
         variable_type='real',
         variable_boundaries=varbound_T,
          algorithm_parameters=algorithm_param)

modelT.run()

Логарифмическая вероятность 8150, уже лучше, чем нормальное распределение, теперь мы попробуем распределение Коши:

def Cauchy_likelihood(params):
gamma = params[1]
loc = params[0]
data = log_rets
log_like = cauchy.logpdf(np.array(data), loc = loc, scale = gamma).sum()
return -log_like

varbound_cauchy= np.array([[-1,1],[0,10]])
algorithm_param = {'max_num_iteration': 1000,\
'population_size':1000,\
'mutation_probability':0.2,\
'elit_ratio': 0.05,\
'crossover_probability': 0.5,\
'parents_portion': 0.3,\
'crossover_type':'uniform',\
'max_iteration_without_improv':None}

modelcauchy=ga(function=Cauchy_likelihood,
dimension=2,
variable_type='real',
variable_boundaries=varbound_cauchy,
algorithm_parameters=algorithm_param)

modelcauchy.run()

Теперь стабильный выпуск:

def Levy_stable_log_likelihood(params):
alpha, beta, scale, loc = params[0], params[1], params[2], params[3]
data = log_rets
Log_likelihood = np.sum(stats.levy_stable.logpdf(np.array(data), alpha, beta,loc = loc, scale = scale))
return -Log_likelihood

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

from pyswarm import pso

lb = [0,-1,0,-0.1]
ub = [2,1,0.1,0.1]

xopt, fopt = pso(Levy_stable_log_likelihood, lb, ub,
swarmsize = 50,
maxiter = 100,
omega = 2,
phip = 2,
phig = 1,
minfunc = 1e-1,
debug = True
)

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

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

import GPyOpt

def Levy_stable_log_likelihood(params):
alpha, beta, scale, loc = params[:,0], params[:,1], params[:,2], params[:,3]
data = log_rets
Log_likelihood = np.sum(levy_stable.logpdf(np.array(data), alpha, beta,loc = loc, scale = scale))
return -Log_likelihood


space =[{'name': 'alpha', 'type': 'continuous', 'domain': (1.51938962,2)},
{'name': 'beta', 'type': 'continuous', 'domain': (-0.8,-0.4)},
{'name': 'scale', 'type': 'continuous', 'domain': (1e-6,0.1)},
{'name': 'loc', 'type': 'continuous', 'domain': (-0.1,0.1)}
]

bo = GPyOpt.methods.BayesianOptimization(f = Levy_stable_log_likelihood,
domain = space,
acquisition_type='EI',
exact_feval = True,
evaluator_type = 'random',
normalize_Y = True,
verbosity = True,
)
bo.run_optimization(max_iter=100, max_time=1000, eps=1e-8, verbosity=True)
cost = bo.fx_opt
points = bo.x_opt
print('Cost function: {} at points {}'.format(cost, points))

Теперь мы можем провести небольшое сравнение максимального логарифмического правдоподобия, чтобы определить наиболее подходящее распределение:

Здесь мы ясно видим, что распределение Стьюдента-t принимает приз как наиболее подходящее распределение для доходности этой акции, с тесным продолжением в стабильном распределении.

Заметка: Этот результат не подлежит обобщению, можно намеренно выборить возвраты таким образом, чтобы результаты отличались.

Чтобы закончить эту статью, я поместил здесь небольшой визуальный элемент для сравнения различных дистрибутивов; распределение Коши было опущено из-за того, что оно было слишком жирным (оно предсказывало ежедневную доходность выше 200% и ниже 100%).

Часть 3: Куда вы идете с этим?

Ну, я не могу сказать вам нигде конкретно, если это то, что вы хотите, но я могу дать вам диапазон, используя моделирование Монте-Карло.

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

Засучив рукава

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

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

Здесь мы сразу же и наиболее зримо наблюдаем недостатки гауссовского соответствия: экстремальные доходности сильно недооцениваются, для дальнейшего сравнения мы сравниваем функции выживания каждого распределения, оцененные с положительной доходностью 20%:

Нормальное распределение недооценивает экстремальные явления на несколько порядков, например, оно недооценивает вероятность доходности, превышающей или равной 20% в любой данный день, примерно в 10 раз в 17-й степени.

Прежде чем мы приготовим

Как я уже говорил в начале этой статьи, я не могу сказать вам конкретно, куда это нас приведет, но что это значит? Это означает, что моделирование методом Монте-Карло не является предсказанием, оно не делает точных заявлений о будущем. Они предназначены только для имитации определенного стохастического процесса, а это означает, что неопределенность присуща вашим результатам. Более того, результаты моделирования сильно зависят от предположений, которые вы делаете о моделируемом процессе. Например, здесь мы неявно предположили, что доходность равна:

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

Это замечательный день, чтобы быть спагетти.

Во-первых, нам нужно решить две вещи:

  1. Сколько сценариев мы хотим изучить?
  2. Через сколько дней в будущем будут смоделированы наши сценарии?

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

Для начала, предположим, что мы хотим увидеть, где мы могли бы оказаться в течение следующих 2 недель, если бы мы купили ETEL (биржевой тикер, с которым мы работаем с начала серии) по какой-то цене сегодня, мы будем сэмплировать с помощью rvs Scipy.stats (random variables):

# Distribution parameters
v = 2.91162520
scale = 1.21872294e-02
loc = -8.13209384e-05
# importing the datetime library to generate dates for the future

from datetime import datetime, timedelta
from scipy import stats

# the simulations' starting price
initial_price = ETEL_clean['Price'].iloc[-1]

# number of simulation days
N_steps = 14
# number of scenarios in simulation
N_sims = 500000
# here we sample the distribution for 14x50000 future returns
rets_t= stats.t.rvs(df= v,loc = loc, scale =scale, size = (N_steps, N_sims))
# here we calculate the price at each timestep from the generated log returns
sim = (np.exp(rets_t.cumsum(axis = 1)))
# we save the results in a dataframe and put them on equal footing
results = pd.DataFrame(np.array(sim).reshape(N_sims, N_steps).T)
results = results.divide(results.iloc[0])*initial_price

time = pd.date_range(start = ETEL_clean.index[-1],
end = ETEL_clean.index[-1] + timedelta(days = N_steps - 1))
# lets look at the final price
final_price_t = results.iloc[-1]
plt.hist(final_price_t[final_price_t < 5*initial_price], density = True, bins = 10000)
plt.title(f'Future Price distribution after {N_steps} days over {N_sims} scenarios')
plt.xlabel('Price/ EGP')
plt.ylabel('Probability Density')

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

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

Давайте посмотрим на распределение цен на целый год вперед (из-за ограничений вычислительной мощности на моем скромном ноутбуке количество сценариев будет ограничено 100000):

Хорошие спагетти, я прав? Здесь мы сразу видим, что плотность вероятности будущей цены асимметрична, это потому, что она ограничена так, что цена неотрицательна, поэтому она может упасть только настолько, насколько сильно. Мы также сразу видим, что кривая вероятности стала грубой и неровной, это потому, что мы уменьшили количество сценариев, поэтому симуляция не совсем сошлась с требуемой нам формой. Мы можем увеличить количество сценариев до 500000, но я не смогу показать историю цен, потому что для отображения требуется 30 минут:

Здесь мы имеем гораздо более плавную кривую, ближе к ожидаемому распределению логарифмического Стьюдента.

Что это говорит нам?

Прежде чем мы ответим на этот вопрос, мы должны сначала ответить на его отрицательный аналог: что это не говорит нам?

  1. Он не говорит нам, есть ли на самом деле низкое или высокое изменение, есть ли низкая или высокая вероятность того, что цена пойдет куда-либо в течение определенного периода. Например, если есть недавние новости, которые негативно влияют на ETEL, очевидно, что более вероятно, что акции упадут (внезапное изменение), чем предполагает распределение.
  2. Он не говорит нам о самой низкой возможной цене, другими словами, только потому, что 14-дневная плотность вероятности показывает, что 18 EGP (например) является самой низкой ценой на гистограмме, не означает, что акция не может опуститься ниже этого, она определенно может, даже при самых обыденных рыночных условиях)
  3. Это повторение первого пункта: то, что гистограмма достигает пика по определенной цене, не означает, что это то, чего вы должны ожидать. Существует бесконечно много путей для акций, и они не обязаны подчиняться распределению вероятностей. Доходность акций имеет много статистических свойств, которые НЕ включены в эту симуляцию.
  4. Мы не просто меняем цвет платья с нормального на какое-то распределение жирного хвоста, прекрасно знаем, что эта симуляция не говорит вам ни капли о том, что возможно в будущем. Распределение жира в хвосте вряд ли поддается математическому моделированию.

Вот что он говорит вам:

  1. Не стоит недооценивать вероятность экстремальных явлений, вероятность их возникновения на порядки выше, чем вы можете ожидать.
  2. Это численное приближение для функции плотности вероятности будущих цен. Однако, как указывалось ранее, жирохвостые распределения редко ведут себя достаточно хорошо, чтобы быть значимыми (особенно в финансах, как будет показано в пункте 3)

3. Не стоит недооценивать силу экстремальных событий или даже средних событий. Акция, которая увеличивается на 1% ежедневно в течение целого года, вырастет примерно в 12,6 раза по сравнению с ее стоимостью, это звучит красиво, и все, кроме одного или двух плохих дней с событием «10 сигм» (глядя на вас, LTCM) или даже событием 4 сигмы, может стереть все это. Почти не имеет значения, насколько маловероятны такие события, потому что ожидаемые потери с учетом их наступления слишком катастрофичны.

4. Хотя это моделирование ни в коем случае не является прогностическим, некоторые специалисты используют его для оценки показателей риска, таких как VaR (Value-at-Risk) и CVaR (Conditional Value-at-Risk) за определенный период времени.

Часть 4: Дойная корова Шредингера

В предыдущей статье из серии Monte Carlo Stock Simulation мы успешно провели моделирование по методу Монте-Карло биржевого тикера на египетском фондовом рынке ETEL. В этой статье мы смоделируем портфель акций, состоящий из 6 акций, а затем перейдем к расчету показателей риска на основе этого моделирования.

Краткое резюме

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

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

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

В этой статье мы будем моделировать портфель акций, состоящий из 6 акций, однако представленная здесь процедура может быть распространена на n-активы, которые не обязательно должны быть в одном классе.

В чем отличие?

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

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

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

На большую работу

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

Наш набор данных состоит из 6 историй акций за последние 10 лет:

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

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

from scipy.stats import multivariate_t

def multivariate_t_likelihood(params):
df = params[0]
data = log_returns(m_equal)
size = len(data.columns)
loc = params[1]*data.mean() + params[2]*np.ones(size)
shape = data.cov()*params[3]**2 + data.cov()*params[4] + params[5]*np.eye(size) + params[6]*data.cov()*data.cov()
if np.linalg.det(shape) < 0:
shape = data.cov()*params[3]
log_like = multivariate_t.logpdf(data, loc = loc, shape = shape, df = df)
return -log_like.sum()

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

from geneticalgorithm import geneticalgorithm as ga

varbounds = np.array([
[0,5],
[-1,1],
[-1,1],
[0,1],
[0,1],
[0,1],
[0,1]
])

algorithm_param = {'max_num_iteration': 1000,\
'population_size':100,\
'mutation_probability':0.2,\
'elit_ratio': 0.05,\
'crossover_probability': 0.5,\
'parents_portion': 0.3,\
'crossover_type':'uniform',\
'max_iteration_without_improv':100}
model=ga(function=multivariate_t_likelihood,
dimension=6,
variable_type='real',
variable_boundaries=varbounds,
algorithm_parameters=algorithm_param)

model.run()

Это дает нам решение:

Если мы скептически относимся к этим результатам, мы можем попробовать другой алгоритм, такой как двойной отжиг от Scipy, который дает нам следующие результаты:

Что кажется достаточно близким, однако важно знать, что можно сделать больше, если мы выясним лучшее приближение масштаба и местоположения распределения Стьюдента-t.

Дойная корова Шредингера

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

Теперь мы можем смоделировать будущее каждой акции для произвольного горизонта произвольное количество раз:

Здесь мы видим, что, если присмотреться, существует некоторая корреляция между движениями цен на акции. Однако их корреляция не совсем та корреляция, которую мы наблюдаем в реальной жизни, которая сильно зависит от секторов, в которых находятся эти акции, и величины событий, приводящих к этим изменениям, т.е. зависимость от хвоста. Линейная корреляция — не самое лучшее приближение, но на этот раз мы с ней справимся.

Но нас интересует не это, нас интересует эффективность портфеля в целом, например, мы создаем портфель из этих акций, с произвольным распределением. Затем мы можем смоделировать это портфолио сколько раз, а затем изучить, как изменяются полученные свойства:

N_sims = 10000
N_steps = 252
# here we decide how much of each stock comprises our portfolio
# we choose an equal parts allocation
weights = np.ones(len(data.columns))
weights /= weights.sum()

stocks_returns = multivariate_t.rvs(loc = loc, shape = shape, df = df, size = (N_sims, N_steps))
time = pd.date_range(
start = m_equal.index[-1],
end = m_equal.index[-1] + timedelta(days = N_steps - 1),
freq = '1D'
)
# here we calculate the overall returns
portfolio_returns = stocks_returns*(weights)
portfolio_returns = portfolio_returns.sum(axis = 2)
portfolio_returns = pd.DataFrame(portfolio_returns.T, index = time)
# here we apply the exponent to get portfolio value from the log returns
portfolio_value = portfolio_returns.cumsum().apply(np.exp)
# then we put all the simulations on equal footing
portfolio_value = portfolio_value.divide(portfolio_value.iloc[0])

Давайте посмотрим на итоговое распределение значений:

Здесь мы видим классическое распределение log-t, теперь пришло время рассчитать различные метрики риска.

Познакомьтесь с VaRs

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

В любом случае, давайте перейдем к делу. Аббревиатура VaR расшифровывается как Value-at-Risk, ценность, которую ваши активы имеют при определенном уровне риска. Тем не менее, существует множество VaR, включая, но не ограничиваясь:

  1. VaR (Value-at-Risk) определяется как стоимость ваших активов с определенной вероятностью (обычно от 5% до 1%) наихудших сценариев.
  2. CVaR (условная стоимость под риском) определяется как среднее значение стоимости вашего актива, учитывая, что значение ниже VaR. Эта метрика также требует % шанса в качестве входных данных. Говоря простым языком, если вам настолько не повезло, что с вами случаются только худшие 5% случаев, то CVaR будет вашим средним богатством.
  3. (Даже мое понимание этого шаткое) EVaR (энтропийная ценность в риске) может быть определена как верхняя граница для VaR и CVaR, полученная из неравенства Чернова. EVaR зависит от уровня неприятия риска, параметра, который показывает, насколько инвестор не склонен к риску.

Мы можем легко рассчитать первые два, VaR и CVaR, но для EVaR мы прибегнем к Riskfolio, пакету количественных финансов Python. Для первоначальных инвестиций в размере 2500 EGP показатели косвенного риска рассчитываются следующим образом:

N_sims = 100000
N_steps = 252
initial_investment = 2500
weights = np.ones(len(data.columns))
weights /= weights.sum()
stocks_returns = multivariate_t.rvs(loc = loc, shape = shape, df = df, size = (N_sims, N_steps))
time = pd.date_range(
start = m_equal.index[-1],
end = m_equal.index[-1] + timedelta(days = N_steps - 1),
freq = '1D'
)
portfolio_returns = stocks_returns*(weights)
portfolio_returns = portfolio_returns.sum(axis = 2)
portfolio_returns = pd.DataFrame(portfolio_returns.T, index = time)
portfolio_value = portfolio_returns.cumsum().apply(np.exp)
portfolio_value = portfolio_value.divide(portfolio_value.iloc[0]/initial_investment)
final_value = portfolio_value.iloc[-1]

alpha = 5
def VaR(value, alpha = 5):
return np.percentile(value, alpha)
ValueatRisk = VaR(final_value,alpha)

def CVaR(value, alpha = 5):
VaR = np.percentile(value, alpha)
return np.mean(value[value < VaR])

ConditionalValueatRisk = CVaR(final_value,alpha)

import riskfolio as rp

EntropicValueatRisk, Risk_Aversion = rp.RiskFunctions.EVaR_Hist(final_value,alpha/100)

result = f"Risk Metrics calculated at {N_steps} days\nAt significance 5 % \nValue at Risk: {ValueatRisk:.2f}\nConditional Value at Risk: {ConditionalValueatRisk:.2f}\nEntropic Value at Risk: {-EntropicValueatRisk:.2f}, at risk aversion of {Risk_Aversion:.2f}"
print(result)

Это говорит нам о том, насколько плохими могут быть дела в следующем году, в рамках этой модели мы можем сказать, что в нижних (худших) 5% сценариев портфель упадет примерно до 1770 EGP, если произойдут только худшие 5% случаев, то мы можем ожидать, что портфель упадет до 1615 EGP. Если какая-либо из этих оценок неверна, при максимальной энтропии она упадет до 1450 EGP.

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

Оптимизация параметров и анализ чувствительности для оптимизации качества процесса в Python

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

Рассмотрим производственный процесс, в котором производится критически важный компонент, используемый в промышленном применении. Качество этого компонента имеет жизненно важное значение, так как любые дефекты или отклонения могут привести к сбоям в работе и значительным финансовым потерям. Чтобы обеспечить надежный и надежный процесс, организация стремится оптимизировать качество своего процесса. В настоящее время индекс производительности процесса (Ppk) составляет 0,92, а коэффициент несоответствия спецификациям составляет 0,6%. Для повышения качества процесса крайне важно определить ключевые факторы, влияющие на цель, и контролировать вариабельность входных данных.

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

Ниже приведен регрессионный анализ данных процесса. Y — выходные данные процесса (отклик), а X1, X2, X3 и X4 — входные данные процесса (предикторы) процесса.

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

Тепловая карта корреляции показывает, что корреляция между X2 и X4 составляет -0,97, что указывает на сильную отрицательную линейную зависимость между предикторами. Сильная отрицательная корреляция между двумя предикторами предполагает, что они предоставляют аналогичную информацию модели. В регрессионном анализе эта высокая корреляция между предикторами известна как мультиколлинеарность. Это может повлиять на стабильность коэффициентов регрессии и затруднить интерпретацию индивидуальных эффектов каждого предиктора. Удаление одного из предикторов может помочь уменьшить мультиколлинеарность и упростить интерпретацию модели. Решение о том, какой предиктор удалить, может быть основано на пошаговой регрессии для определения наиболее подходящего подмножества предикторов.

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

Метод обратного исключения исключил предикторы X3 и X4 из модели, в конечном итоге выбрав предикторы X1 и X2 для модели. Это снижает мультиколлинеарность и упрощает интерпретацию модели. Для проверки мультиколлинеарности в регрессионной модели в качестве меры коллинеарности между предикторными переменными можно использовать дисперсионный инфляционный коэффициент (VIF). VIF количественно определяет, насколько завышена дисперсия расчетных коэффициентов регрессии из-за мультиколлинеарности. Значение VIF больше 1 указывает на наличие мультиколлинеарности, а более высокие значения указывают на более сильный эффект коллинеарности.

Значения VIF (1,055) указывают на то, что мультиколлинеарность отсутствует в модели.

Сгенерируйте уравнение регрессии, чтобы оценить взаимосвязь между ответом Y и предикторами X1 и X2.

Коэффициент для X1 в уравнении равен 1,47, что указывает на то, что Y увеличится на 1,47 единицы для увеличения X1 на одну единицу, при условии, что X2 остается постоянным. Это говорит о том, что X1 положительно влияет на значение Y.

Точно так же коэффициент для X2, равный 0,66, предполагает, что при увеличении X2 на одну единицу Y увеличится на 0,66 единицы, предполагая, что X1 остается постоянным. Это указывает на положительную связь между X2 и значением Y.

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

Вот интерпретация ключевых статистических данных:

1. R-квадрат и скорректированный R-квадрат: Значение R-квадрата равно 0,979, что означает, что 97,9% изменчивости Y объясняется X1 и X2. Скорректированный R-квадрат корректирует количество предикторов в модели. Скорректированное значение R-квадрата равно 0,979, что указывает на то, что регрессионная модель хорошо соответствует данным.

2. F-статистика и Prob (F-статистика): F-статистика проверяет общую значимость регрессионной модели. Большая F-статистика (229,5) с небольшим p-значением (4,41e-09) указывает на то, что модель является статистически значимой, предполагая, что по крайней мере один из предикторов значимо связан с ответом.

3. Коэффициенты: Коэффициент const (52.5773) представляет собой пересечение или математическое ожидание Y, когда и X1, и X2 равны нулю. Коэффициенты для X1 (1,4683) и X2 (0,6623) показывают, как изменение на одну единицу каждого предиктора связано с изменением зависимой переменной.

4. P-значения: P-значения, связанные с каждым коэффициентом, проверяют нулевую гипотезу о том, что коэффициент равен нулю. При этом все коэффициенты имеют крайне малые p-значения (близкие к нулю), указывающие на то, что X1 и X2 значимо связаны с Y.

5. Омнибус: Омнибусный тест — это проверка предположения о нормальности остатков. В этом случае статистика Omnibus равна 1,509, а связанное с ней p-значение равно 0,470. Поскольку p-значение больше 0,05, нет никаких существенных доказательств отклонения от нормы в остатках. Это указывает на то, что предположение о нормальности является разумным для модели.

6. Жарк-Бера: Тест Жарка-Бера — еще одна проверка предположения о нормальности остатков. Он исследует как асимметрию, так и эксцесс, чтобы оценить, следуют ли остатки нормальному распределению. В этом случае статистика Жарка-Бера равна 1,104, а связанное с ней p-значение равно 0,576. Как и в случае с омнибусным тестом, p-значение превышает 0,05, что указывает на отсутствие существенных признаков отклонения от нормы.

7. Дурбин-Уотсон: Статистика Дурбина-Уотсона используется для обнаружения наличия автокорреляции (корреляции остатков во времени) в модели. Значение Дурбина-Уотсона, равное 1,922, указывает на положительную автокорреляцию, что означает, что остатки положительно коррелируют с их запаздывающими значениями. Значение от 1 до 2 предполагает мягкую положительную автокорреляцию.

Хотя статистические тесты, такие как Омнибус, Жарк-Бера и Дурбин-Уотсон, могут дать представление о нормальности остатков и автокорреляции остатков, они могут не дать подробного представления о природе и величине отклонений. Чтобы получить более полное понимание, гистограммы, графики Q-Q и диаграммы рассеяния будут использоваться в сочетании со статистическими тестами для тщательной оценки предположений.

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

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

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

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

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

https://medium.com/@chenchungwai/a-comprehensive-guide-to-regression-model-validation-with-visual-tools-and-statistical-tests-in-f86d4ff6efbb

В целом, модель, по-видимому, хорошо подходит для данных. Результаты показывают, что X1 и X2 имеют сильную положительную связь с Y, на что указывают положительные коэффициенты и их статистическая значимость. Остатки в модели распределены приблизительно нормально без значительных отклонений от нормальности.

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

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

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

Ниже приведены шаги для выполнения моделирования методом Монте-Карло:

  1. Определите распределение вероятностей входных переменных:

Используйте исторические данные для определения среднего значения и стандартного отклонения входных переменных X1 и X2:

Среднее значение X1: 7,5
Стандартное отклонение X1: 0,6
Среднее значение X2: 48
Стандартное отклонение X2: 1,6

2. Определите количество симуляций для выполнения:

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

В этом случае мы выполним 1 000 000 симуляций.

Во время каждого моделирования случайные выборки для X1 и X2 генерируются на основе их соответствующих распределений вероятностей с использованием заданных средних значений и стандартных отклонений. Выходная переменная (Y) вычисляется по уравнению регрессии Y = 52,58 + 1,47 * X1 + 0,66 * X2. Значения выходных переменных (Y), полученные в результате каждого моделирования, собираются для анализа распределения значений выходных переменных. Сравнивая смоделированные значения Y с заданным целевым значением (95,4) и нижним и верхним пределами спецификации (92,4 и 98,4), мы можем оценить индекс производительности процесса (Ppk) и коэффициент отклонения от спецификации. Этот анализ помогает нам определить вероятность соответствия спецификациям и выявить любые потенциальные проблемы с качеством процесса.

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

Это скрипт Python для моделирования Монте-Карло.import numpy as np

# Define the regression equation coefficients
intercept = 52.58
coef_x1 = 1.468
coef_x2 = 0.6623

# Define the target value, lower spec, and upper spec for y
target_y = 95.4
lower_spec = 91.6
upper_spec = 99.2

# Define the mean and standard deviation for X1 and X2
mean_x1 = 7.5
std_x1 = 0.6
mean_x2 = 48
std_x2 = 1.3

# Set the number of iterations for the Monte Carlo simulation
num_iterations = 1000000

# Initialize array to store the simulated y values
simulated_y = np.zeros(num_iterations)

# Perform the Monte Carlo simulation
for i in range(num_iterations):
# Generate random samples for X1 and X2
x1 = np.random.normal(mean_x1, std_x1)
x2 = np.random.normal(mean_x2, std_x2)

# Calculate the predicted value of y using the regression equation
y_predicted = intercept + coef_x1 * x1 + coef_x2 * x2

# Store the predicted value of y
simulated_y[i] = y_predicted

# Calculate the mean and standard deviation of the simulated y values
mean_simulated_y = np.mean(simulated_y)
std_simulated_y = np.std(simulated_y)

# Calculate the count and percentage of out-of-specification values
out_of_spec_count = np.sum((simulated_y < lower_spec) | (simulated_y > upper_spec))
out_of_spec_percentage = out_of_spec_count / num_iterations * 100

# Calculate the process capability (Cp)
cp = (upper_spec — lower_spec) / (6 * np.std(simulated_y))

# Calculate the process centering (Cpk)
cpk = min((upper_spec — mean_simulated_y) / (3 * np.std(simulated_y)),
(mean_simulated_y — lower_spec) / (3 * np.std(simulated_y)))

# Calculate the Process Performance Index (Ppk)
ppk = min(cp, cpk)

# Print the simulation results
print(«Mean of simulated y values:», mean_simulated_y)
print(«Standard deviation of simulated y values:», std_simulated_y)
print(«Out-of-Specification Count:», out_of_spec_count)
print(«Out-of-Specification Percentage:», out_of_spec_percentage)
print(«Process Capability (Cp):», cp)
print(«Process Centering (Cpk):», cpk)
print(«Process Performance Index (Ppk):», ppk)

Это результаты, полученные с помощью моделирования методом Монте-Карло

Среднее значение составляет 95,38, что указывает на то, что среднее прогнозируемое значение Y близко к целевому значению 95,4.
Стандартное отклонение составляет 1,378, что говорит о некоторой изменчивости прогнозируемых значений Y.
Количество значений, не соответствующих спецификации, составляет 58453, что указывает на то, что из 1 миллиона симуляций смоделировано 58453 значений Y, которые выходят за указанный диапазон нижнего и верхнего пределов спецификации.
Процент отклонений от спецификации составляет 0,58%, что дает представление о доле смоделированных значений, которые не соответствуют указанным требованиям к качеству.
Значение Cp составляет 0,919, что говорит о том, что процесс имеет некоторую потенциальную возможность соответствовать спецификациям.
Значение Cpk составляет 0,914, что указывает на то, что процесс немного смещен в сторону одного из пределов спецификации.
Значение Ppk составляет 0,914, что означает, что общая производительность процесса составляет 0,914.

Результаты моделирования, полученные с помощью Monte Carlo Simulation, точно соответствуют текущим характеристикам процесса. Текущий индекс производительности процесса (Ppk) соответствует смоделированному значению PPK, что указывает на то, что моделирование точно отражает поведение процесса. Однако, поскольку наша цель состоит в том, чтобы достичь значения Ppk 1,67 или выше, мы признаем необходимость улучшения процесса. Чтобы определить конкретные входные переменные, которые достигают этой цели, мы обратимся к оптимизации параметров. Оптимизация параметров — это метод, используемый для нахождения оптимальных значений входных переменных процесса, которые приводят к желаемым результатам или целям. Используя методы оптимизации, мы можем определить оптимальную комбинацию входных переменных, которая приведет к желаемому значению PPK.

Ниже приведены шаги по оптимизации параметров.

  1. Определите диапазон X1 и X2: установите допустимый диапазон или границы, в пределах которых могут изменяться значения X1 и X2.

Диапазон X1 составляет от 7,3 до 7,7
Диапазон X2 составляет от 47,8 до 48,2

2. Определите целевую функцию: Целевой функцией является максимизация индекса производительности процесса (Ppk).

3. Выполните моделирование методом Монте-Карло: В цикле с указанным числом итераций, установленным на 1 000 000, генерируются случайные выборки X1 и X2 в определенных диапазонах. Прогнозируемое значение y вычисляется с помощью уравнения регрессии. Вычисляется значение Cpk для моделируемого среднего значения y. Если вычисленное значение Cpk превышает текущее максимальное значение Cpc, максимальное значение Cpk и соответствующие значения X1 и X2 обновляются.

Выполняя оптимизацию параметров, мы можем определить оптимальные значения X1 и X2, которые максимизируют индекс производительности процесса (Ppk)

Это скрипт Python для оптимизации параметров.# Define the range of X1 and X2
x1_range = [7.3, 7.7]
x2_range = [47.8, 48.2]

# Set the number of iterations for the Monte Carlo simulation
num_iterations = 1000000

def calculate_cpk(mean_y, std_y):
cp = (upper_spec — lower_spec) / (6 * std_y)
cpk = min((upper_spec — mean_y) / (3 * std_y),
(mean_y — lower_spec) / (3 * std_y))
return cpk

# Initialize variables to store the highest Cpk and its corresponding X1 and X2
highest_cpk = 0.0
best_x1 = 0.0
best_x2 = 0.0

# Perform Monte Carlo Simulation
for i in range(num_iterations):
# Generate random samples for X1 and X2 within the specified range
x1 = np.random.uniform(*x1_range)
x2 = np.random.uniform(*x2_range)

# Calculate the predicted value of y using the regression equation
y_predicted = intercept + coef_x1 * x1 + coef_x2 * x2

# Calculate the standard deviation of y assuming a normal distribution
std_y = (upper_spec — lower_spec) / (6 * np.sqrt(3))

# Calculate the Cpk value for the simulated mean value of y
cpk = calculate_cpk(y_predicted, std_y)

# Update the highest Cpk and its corresponding X1 and X2 if necessary
if cpk > highest_cpk:
highest_cpk = cpk
optimized_x1 = x1
optimized_x2 = x2

# Calculate additional values based on the highest Cpk, optimized X1, and optimized X2
mean_simulated_y = intercept + coef_x1 * optimized_x1 + coef_x2 * optimized_x2
std_simulated_y = (upper_spec — lower_spec) / (6 * np.sqrt(3))
out_of_spec_count = np.sum((simulated_y < lower_spec) | (simulated_y > upper_spec))
out_of_spec_percentage = out_of_spec_count / num_iterations * 100
cp = (upper_spec — lower_spec) / (6 * std_simulated_y)
ppk = min(cp, highest_cpk)

# # Print the highest Cpk value and its corresponding X1 and X2
print(«Highest Cpk value:», highest_cpk)
print(«Optimized X1 value:», optimized_x1)
print(«Optimized X2 value:», optimized_x2)
# Print the optimized results
print(«Optimized Results:»)
print(«Mean of simulated y values:», mean_simulated_y)
print(«Standard deviation of simulated y values:», std_simulated_y)
print(«Out-of-Specification Count:», out_of_spec_count)
print(«Out-of-Specification Percentage:», out_of_spec_percentage)
print(«Process Capability (Cp):», cp)
Print(«Process Capability Index (Cpk):», highest_cpk)
print(«Process Performance Index (Ppk):», ppk)

Это результаты, полученные при оптимизации параметров

Процесс оптимизации параметров привел к значительному улучшению качества процесса, о чем свидетельствуют следующие результаты:

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

Оптимизированное значение X1 равно 7,46, а оптимизированное значение X2 — 48,13. Эти значения представляют собой оптимальные настройки входных переменных, которые максимизируют производительность процесса.

Среднее значение смоделированных значений y составляет 95,40, что соответствует целевому значению 95,4. Это указывает на то, что процесс был доработан для достижения желаемого результата.

Стандартное отклонение моделируемых значений y составляет 0,73, что свидетельствует о значительном снижении вариабельности выходных данных.

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

Процент несоответствия спецификациям составляет 0,0%.

Значение возможности процесса (Cp) теперь равно 1,73, процесс центрируется, когда Cp совпадает с Cpk.

Значение индекса производительности процесса (Ppk) соответствует самому высокому значению Cpk 1,73, что отражает общую производительность процесса и его способность постоянно соответствовать спецификациям.

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

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

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

Основываясь на выводах, полученных в результате оптимизации параметров, мы теперь обратимся к анализу чувствительности, чтобы оценить влияние изменения стандартных отклонений входных переменных X1 и X2 на процент значений, отклоняющихся от спецификации, для выходной переменной Y.

Вот шаги для выполнения анализа чувствительности:

1. Определите диапазон процентных изменений для стандартных отклонений. Этот диапазон определяет степень, в которой будут скорректированы стандартные отклонения.

Диапазон процентных изменений для стандартных отклонений составляет от -20% до 20%.

2. Выполните моделирование методом Монте-Карло, сгенерировав случайные выборки для X1 и X2, корректируя их стандартные отклонения на основе определенных диапазонов.

3. Рассчитайте процент значений, отклоняющихся от спецификации, как для X1, так и для X2. Это позволяет количественно оценить влияние изменения стандартных отклонений на значения, не соответствующие спецификации.

4. Постройте диаграмму рассеяния, чтобы визуализировать результаты анализа чувствительности.

Это скрипт Python для анализа чувствительности.# Define the range of percentage changes for standard deviations
percentage_changes = np.linspace(-0.20, 0.20, num=9)

# Initialize lists to store the out-of-specification percentages for X1 and X2
out_of_spec_percentage_x1 = []
out_of_spec_percentage_x2 = []

# Perform the Monte Carlo simulation for different standard deviations of X1 while holding X2 constant
for percentage_change in percentage_changes:
# Calculate the adjusted standard deviation for X1
std_x1_adjusted = std_x1 * (1 + percentage_change)

# Initialize counter for out-of-specification instances
out_of_spec_count_x1 = 0

# Perform the Monte Carlo simulation
for i in range(num_iterations):
# Generate random sample for X1 with fixed values
x1 = np.random.normal(mean_x1, std_x1_adjusted)

# Calculate the predicted value of y using the regression equation
y_predicted = intercept + coef_x1 * x1 + coef_x2 * mean_x2

# Check if the predicted value of y is out of specification for X1
if y_predicted < lower_spec or y_predicted > upper_spec:
out_of_spec_count_x1 += 1

# Calculate the percentage out-of-specification for X1
out_of_spec_percentage_x1.append((out_of_spec_count_x1 / num_iterations) * 100)

# Perform the Monte Carlo simulation for different standard deviations of X2 while holding X1 constant
for percentage_change in percentage_changes:
# Calculate the adjusted standard deviation for X2
std_x2_adjusted = std_x2 * (1 + percentage_change)

# Initialize counter for out-of-specification instances
out_of_spec_count_x2 = 0

# Perform the Monte Carlo simulation
for i in range(num_iterations):
# Generate random sample for X2 with fixed values
x2 = np.random.normal(mean_x2, std_x2_adjusted)

# Calculate the predicted value of y using the regression equation
y_predicted = intercept + coef_x1 * mean_x1 + coef_x2 * x2

# Check if the predicted value of y is out of specification for X2
if y_predicted < lower_spec or y_predicted > upper_spec:
out_of_spec_count_x2 += 1

# Calculate the percentage out-of-specification for X2
out_of_spec_percentage_x2.append((out_of_spec_count_x2 / num_iterations) * 100)

# Plot the scatterplot for both X1 and X2 simulations
plt.scatter(percentage_changes, out_of_spec_percentage_x1, color=’red’, label=’X1′)
plt.plot(percentage_changes, out_of_spec_percentage_x1, color=’red’, linestyle=’-‘, linewidth=2)
plt.scatter(percentage_changes, out_of_spec_percentage_x2, color=’blue’, label=’X2′)
plt.plot(percentage_changes, out_of_spec_percentage_x2, color=’blue’, linestyle=’-‘, linewidth=2)

plt.xlabel(«Change in Percentage of Standard Deviation»)
plt.ylabel(«% Out-of-Specification»)
plt.title(«Simulation Results: % Out-of-Specification vs Change in Standard Deviation»)
plt.legend()
plt.show()

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

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

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

В этой статье я покажу, как можно использовать перспективное моделирование Монте-Карло для прогнозирования будущей стоимости собственного капитала с помощью логнормального распределения.

Фото: Внутренняя стоимость

1. Мотивация

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

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

Упрощенно можно представить себе подход моделирования по методу Монте-Карло как многократное извлечение мячей для гольфа из большой корзины с заменой. Размер и форма корзины зависят от предположения о распределении входных данных (например, нормальное распределение со средним значением 100 и стандартным отклонением 10, по сравнению с равномерным распределением или треугольным распределением), когда некоторые корзины глубже или более симметричны, чем другие, что позволяет вытаскивать определенные шары чаще, чем другие.

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

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

В данной статье показано, как можно использовать прогнозное моделирование по методу Монте-Карло для прогнозирования будущей стоимости собственного капитала с использованием логнормального распределения.

2. Марковское свойство

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

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

Статистические свойства истории цен на акции IBM могут быть полезны для определения характеристик стохастического процесса, за которым следует цена акции (например, ее волатильность). Здесь подчеркивается, что конкретный путь, по которому шли акции в прошлом, не имеет значения.

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

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

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

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

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

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

3. Стохастический процесс непрерывного времени

Рассмотрим переменную, которая следует марковскому стохастическому процессу. Предположим, что его текущее значение равно 10 и что изменение его значения в течение одного года равно φ(0, 1), где φ(μ, σ) обозначает распределение вероятностей, которое обычно распределяется со средним μ и стандартным отклонением σ.

Каково распределение вероятностей изменения значения переменной в течение двух лет? Изменение за два года является суммой двух нормальных распределений, каждое из которых имеет среднее значение, равное нулю, и стандартное отклонение, равное 1,0.

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

Дисперсия распределения вероятностей равна квадрату его стандартного отклонения. Таким образом, дисперсия изменения значения рассматриваемой нами переменной за один год равна 1,0

Таким образом, среднее значение изменения рассматриваемой нами переменной в течение двух лет равно нулю, а дисперсия этого изменения равна 2,0. Таким образом, изменение переменной за два года составляет φ(0, √2).

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

Мы предполагаем, что это одно и то же. Из этого следует, что дисперсия изменения в течение шестимесячного периода должна составлять 0,5. Эквивалентно, стандартное отклонение изменения равно √0,5, так что распределение вероятностей изменения значения переменной в течение шести месяцев равно φ(0, √0,5)

Аналогичный аргумент показывает, что изменение значения переменной в течение трех месяцев равно φ(0, √0,25). В более общем случае, изменение в течение любого периода времени длины T составляет φ(0, √T).

В частности, изменение длины δt за очень короткий период времени составляет φ(0, √δt). Знаки квадратного корня в этих результатах могут показаться странными.

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

Дисперсия изменения переменной в нашем примере равна 1,0 в год, так что дисперсия изменения за два года равна 2,0, а дисперсия изменения за три года равна 3,0. Стандартное отклонение изменения через два и три года составляет √2 и √3 соответственно.

Строго говоря, мы не должны говорить о стандартном отклонении переменной как о 1,0 в год. Это должно быть «1,0 на квадратный корень из лет». Полученные результаты объясняют, почему неопределенность часто называют пропорциональной квадратному корню из времени.

2.1. Винские процессы

Процесс, за которым следует рассматриваемая нами переменная, известен как винеровский процесс. Это особый тип марковского стохастического процесса со средним изменением, равным нулю, и коэффициентом дисперсии 1,0 в год.

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

Выражаясь формально, переменная z следует винеровскому процессу, если она обладает следующими двумя свойствами:

Свойство 1. Изменение δz за небольшой промежуток времени δt равно

где ε — случайный рисунок из стандартизированного нормального распределения, φ(0, 1).

Свойство 2. Значения δz для любых двух различных коротких интервалов времени δt независимы.

Из первого свойства следует, что δz само по себе имеет нормальное распределение с

Второе свойство подразумевает, что z следует марковскому процессу. Рассмотрим увеличение величины z в течение относительно длительного периода времени, T.

Это можно обозначить как z(T) — z(0). Его можно рассматривать как сумму приращений z за N малых интервалов времени длиныδ t, где

Таким образом

где ε(i = 1,2, . . . , N) — случайные рисунки из φ(0, 1). Из второго свойства винеровских процессов εi независимы друг от друга. Из приведенного выше уравнения следует, что z(T) — z(0) нормально распределяется с

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

По истечении пяти лет он обычно распределяется со средним значением 25 и стандартным отклонением √5, или 2,236. Наша неопределенность относительно значения переменной в определенный момент времени в будущем, измеряемая ее стандартным отклонением, возрастает как квадратный корень из того, как далеко мы смотрим вперед.

В обычном исчислении обычно исходят от небольших изменений к пределу по мере того, как малые изменения приближаются к нулю. Таким образом, δy/δx становится dy/dx в пределе, и так далее.

Аналогично можно поступить и при работе со стохастическими процессами. Винеровский процесс является пределом в виде δt → 0 описанного выше процесса для z.

На рисунке ниже показано, что происходит с траекторией, по которой следует z, при приближении к пределу δt → 0. Обратите внимание, что путь довольно «неровный».

Как получается винеровский процесс при δt → 0

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

  1. Ожидаемая длина пути, за которым следует z в любом интервале времени, бесконечна.
  2. Ожидаемое количество раз, когда z равно какому-либо конкретному значению в любом интервале времени, бесконечно.

2.2. Обобщенный винеровский процесс

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

Коэффициент дисперсии, равный 1,0, означает, что дисперсия изменения z во временном интервале длины T равна T. Обобщенный винеровский процесс для переменной x может быть определен в терминах dz следующим образом:

где a и b — константы.

Чтобы понять приведенное выше уравнение, полезно рассмотреть два компонента в правой части по отдельности. Член a dt подразумевает, что x имеет ожидаемую скорость дрейфа a в единицу времени.

Без члена b dz уравнение имеет вид

из чего следует, что

Интегрируя по времени, получаем

где x0 — значение x в нулевой момент времени. За период времени длины T величина x возрастает на величину aT.

Член b dz в правой части уравнения можно рассматривать как добавление шума или изменчивости к пути, по которому следует x. Величина этого шума или изменчивости равна b, умноженной на винеровский процесс.

Винский процесс имеет стандартное отклонение, равное 1,0. Из этого следует, что b, умноженное на винеровский процесс, имеет стандартное отклонение b.

На малом интервале времени δt изменение значения δx задается приведенными выше уравнениями как

где, как и раньше, ε — случайный отбор из стандартизированного нормального распределения. Таким образом, δx имеет нормальное распределение с

Аргументы, аналогичные аргументам, приведенным для винеровского процесса, показывают, что изменение значения x в любом интервале времени T нормально распределяется с

Таким образом, обобщенный винеровский процесс, приведенный в приведенном выше уравнении, имеет ожидаемый коэффициент дрейфа (т.е. средний дрейф в единицу времени) a и коэффициент дисперсии (т.е. дисперсию в единицу времени) . Это проиллюстрировано на рисунке ниже

Обобщенный винеровский процесс: a = 0,3, b = 1,5

Рассмотрим ситуацию, когда денежная позиция компании, измеряемая в тысячах долларов, подчиняется обобщенному процессу Винера с дрейфом 20 в год и коэффициентом отклонения 900 в год.

Изначально денежная позиция равна 50. В конце одного года денежная позиция будет иметь нормальное распределение со средним значением 70 и стандартным отклонением √900 или 30.

По истечении шести месяцев он будет иметь нормальное распределение со средним значением 60 и стандартным отклонением 30√0,5 = 21,21. Наша неопределенность в отношении денежной позиции в какой-то момент в будущем, измеряемая ее стандартным отклонением, возрастает по мере того, как квадратный корень из того, насколько далеко вперед мы смотрим.

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

2.3. Процесс Ито

Можно определить еще один тип стохастического процесса. Это называется процессом Ито. Это обобщенный винеровский процесс, в котором параметры a и b являются функциями величины базовой переменной x и времени t.

Алгебраически процесс Itô может быть записан

Как ожидаемая скорость дрейфа, так и степень дисперсии процесса Itô могут изменяться с течением времени. В малом промежутке времени между t и t + δt переменная изменяется от x до x + δx, где

Эта зависимость включает в себя небольшое приближение. Предполагается, что дрейф и дисперсия x остаются постоянными, равными a(x, t) и b(x, t)² соответственно, в течение интервала времени между t и t + δt.

3. Модель поведения цены акций

Модель поведения цены акций известна как геометрическое броуновское движение. Дискретно-временная версия модели

Переменная δS — это изменение цены акции S за небольшой интервал времени δt, а ε — случайная выборка из стандартизированного нормального распределения (т.е. нормального распределения со средним значением нуля и стандартным отклонением 1,0).

Параметр μ — это ожидаемая норма доходности в единицу времени от акции, а параметр σ — волатильность цены акции. Оба эти параметра считаются постоянными. Левая часть следующего уравнения

– доходность, обеспечиваемая запасом за короткий промежуток времени δt. Член μδt — это математическое ожидание этой доходности, а член σεδt — стохастическая составляющая доходности.

Дисперсия стохастической составляющей (и, следовательно, всей доходности) равна σ²δt. Это согласуется с определением волатильности как данности; То есть, σ такова, что σ√δt является стандартным отклонением доходности за короткий промежуток времени δt.

Приведенные выше уравнения показывают, как dS/S нормально распределяется со средним значением μ δt и стандартным отклонением σδt. Другими словами,

4. Свойство Lognormal

Теперь мы используем лемму Ито для получения процесса, за которым следует ln S, когда S следует следующему процессу

Определять

При постоянном μ и σ является разумной моделью движения цен на акции. Из леммы Ито следует, что процесс, за которым следует функция G из S и t, равен

Заметим, что и S, и G подвержены влиянию одного и того же основного источника неопределенности, dz. Это оказывается очень важным при выводе результатов Блэка-Шоулза.

Потому что

Из приведенного выше уравнения следует, что процесс, которому следует G, равен

Поскольку μ и σ постоянны, это уравнение показывает, что G = ln S подчиняется обобщенному винеровскому процессу. Он имеет постоянную скорость дрейфа μ —σ²/2 и постоянную дисперсию σ².

Таким образом, изменение ln S между нулевым временем и некоторым будущим временем, T, нормально распределяется со средним значением

и дисперсия

Это означает, что

где ST — цена акции в будущий момент времени TS0 — цена акции в нулевой момент времени, а φ(m, s) обозначает нормальное распределение со средним m и стандартным отклонением s .

Следующее уравнение

показывает, что ln ST нормально распределена. Переменная имеет логнормальное распределение, если натуральный логарифм переменной нормально распределен.

Таким образом, модель поведения цены акций, которую мы разработали в этой статье, подразумевает, что цена акции в момент времени T, учитывая ее сегодняшнюю цену, логарифмически распределена. Стандартное отклонение логарифма цены акции составляет σt. Он пропорционален квадратному корню из того, как далеко вперед мы смотрим.

5. Реальный мир и нейтральный к риску мир

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

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

К сожалению, нелегко узнать правильную ставку дисконтирования, чтобы применить ее к ожидаемому выигрышу в реальном мире. Позиция по колл-опциону более рискованна, чем позиция по акциям, и предположим, что в реальном мире ожидаемая доходность по акциям составляет 16%.

В результате, ставка дисконтирования, применяемая к выплате по колл-опциону, превышает 16%. Не зная стоимости опциона, мы не знаем, насколько она должна быть больше 16%.

Риск-нейтральная оценка решает эту проблему. Мы знаем, что в мире с нейтральным риском ожидаемая доходность всех активов (и, следовательно, ставка дисконтирования, используемая для всех ожидаемых выплат) является безрисковой ставкой.

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

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

Для некоторого значения рыночной цены риска мы получаем «реальный мир» и наблюдаемые на практике темпы роста цен ценных бумаг.

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

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

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

Из этого следует, что актуарный подход верен только тогда, когда математическое ожидание от дериватива одинаково как в реальном, так и в нейтральном к риску мире.

Мы упоминали, что когда мы переходим из реального мира в мир нейтральный к риску, волатильность переменных остается прежней, но ожидаемые темпы их роста могут измениться. Например, ожидаемые темпы роста индекса фондового рынка уменьшаются, возможно, на 4% или 5%, когда мы переходим из реального мира в мир нейтральный к риску.

Кроме того, можно обоснованно предположить, что ожидаемый темп роста переменной одинаков как в реальном, так и в нейтральном к риску мире, если переменная имеет нулевой систематический риск (β=0), так что процентные изменения переменной имеют нулевую корреляцию с доходностью фондового рынка. Из этого можно сделать вывод, что актуарный подход к оценке производного инструмента дает правильный ответ, если все базовые переменные имеют нулевой систематический риск.

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

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

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

6. Перспективное моделирование по методу Монте-Карло

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

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

  1. Оцените ожидаемую доходность капитала в реальном мире для горизонта прогнозирования.
  2. Выборка случайного пути для S в реальном мире.
  3. Повторите шаг 2, чтобы получить множество примеров будущей стоимости акций в реальном мире.
  4. Вычислите среднее значение выборки будущей стоимости собственного капитала, чтобы получить оценку ожидаемой будущей стоимости собственного капитала.

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

где dz — это процесс Винера, μ — ожидаемая доходность в реальном мире, а σ — волатильность. Чтобы смоделировать путь, по которому следует S, мы можем разделить горизонт прогноза на N коротких интервалов длины Δt и получим

где S(t) обозначает значение S в момент времени tε — случайная выборка из нормального распределения со средним нулем и стандартным отклонением 1,0. Это позволяет вычислить значение S в момент времени Δ t из начального значения S, значение в момент времени 2Δ t вычислить из значения в момент времени Δ t и т. д.

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

Из леммы Ито процесс, за которым следует ln S, имеет вид

Таким образом, чтобы

или-эквивалентно

Это уравнение используется для построения пути для S.

Работа с ln S, а не с S дает большую точность. Кроме того, если μ и σ постоянны, то

справедливо для всех Т. Из этого следует, что

Это уравнение может быть использовано для прогнозирования будущей стоимости собственного капитала в момент времени T.

7. Пример перспективного моделирования методом Монте-Карло с использованием логнормального распределения

Предположим, что ожидаемая доходность собственного капитала (рассчитанная на основе CAPM) составляет 15% годовых (в реальном мире), а стандартное отклонение доходности (т.е. волатильность) составляет 50% годовых. Приведенная стоимость капитала в момент времени 0 равна 100 долл. Спрогнозируйте будущую стоимость акций через 5 лет и 95% доверительный интервал для будущей стоимости акций через 5 лет.

Это означает, что S0 = 100,00, μ = 0,15, σ = 0,50 и T = 5,00.

Во-первых, давайте загрузим соответствующие библиотеки pythonimport numpy as np
import pandas as pd
from scipy.stats import norm
from numpy.random import randn
from numpy import random as rn
from matplotlib import pyplot as plt
import seaborn as sns
%matplotlib inline

Во-вторых, давайте напишем наши входные переменныеS0 = 100.00
μ = 0.15
σ = 0.50
T = 5

В-третьих, давайте найдем будущую стоимость капитала через 5 лет в детерминированном мире (т.е. мире абсолютной определенности)

ST = S0*np.exp(μ*T)
ST

211.7000016612675

В-четвертых, решим использовать 50 000 путей и независимое от пути моделирование (т.е. N = 1, то есть разделим горизонт прогноза на N коротких интервалов длины Δt)Paths = 50000
N = 1
dt = T/N

В-пятых, давайте найдем будущую стоимость акций через 5 лет в стохастическом мире (т.е. мире неопределенности). Итак, мы начинаем со столбца под названием S0то есть с приведенной стоимости собственного капитала в момент времени 0data = S0*np.ones(Paths)
# initialize list of lists
data1 = {‘S0’:data}
# Create the pandas DataFrame
df = pd.DataFrame(data1)
# print dataframe.
df

Затем мы создадим столбец с именем RAND(), который генерирует случайные числа от 0 до 1.df[‘RAND()’] = rn.rand(Paths,1)
df

После этого мы создадим столбец с именем NORMSINV(), который дает обратную кумулятивную функцию для стандартного нормального распределения. Из этого следует, что NORMSINV(RAND()) дает случайную выборку из стандартного нормального распределения.df[‘NORMSINV()’] = norm.ppf(df[‘RAND()’])
df

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

df[‘ST’] = df[‘S0’]*np.exp((μ-0.5*σ**2)*T+σ*df[‘NORMSINV()’]*np.sqrt(T))
df

В-шестых, давайте проанализируем информацию моделированияV = df[‘ST’]

G = 1629562571
np.random.seed(G)
from datetime import datetime
start_time = datetime.now()
from datetime import datetime
# datetime object containing current date and time
now = datetime.now()
dt_string = now.strftime(«%d/%m/%Y %H:%M:%S»)
def NORMSINV(x):
x = si.norm.ppf(x)
return (x)
time_elapsed1 = datetime.now() — start_time
from scipy import stats
Workbook_Name = «Prospective MC Simulation.ipynb»
Number_of_Simulations = «{:,.0f}».format(Paths)
Number_of_Iterations = «{:,.0f}».format(N)
Number_of_Inputs = «{:,.0f}».format(5)
Number_of_Outputs = 38
Sampling_Type = «Latin Hypercube»
Simulation_Start_Time = dt_string
Simulation_Duration = «{}».format(time_elapsed1)
Random_N_Generator = «Mersenne Twister»
Random_Seed = G
e = [«Workbook Name»,»Number of Simulations»,»Number of Iterations»,»Number of Inputs»,»Number of Outputs»,»Sampling Type»,\
«Simulation Start Time»,»Simulation Duration»,»Random # Generator»,»Random Seed»]
f = [Workbook_Name, Number_of_Simulations, Number_of_Iterations, Number_of_Inputs, Number_of_Outputs, Sampling_Type,\
Simulation_Start_Time, Simulation_Duration, Random_N_Generator, Random_Seed]
Per5 = «{:,.4f}».format(np.percentile(V, 5))
P5 = «{:.0%}».format(0.05)
Per10 = «{:,.4f}».format(np.percentile(V, 10))
P10 = «{:.0%}».format(0.10)
Per15 = «{:,.4f}».format(np.percentile(V, 15))
P15 = «{:.0%}».format(0.15)
Per20 = «{:,.4f}».format(np.percentile(V, 20))
P20 = «{:.0%}».format(0.20)
Per25 = «{:,.4f}».format(np.percentile(V, 25))
P25 = «{:.0%}».format(0.25)
Per30 = «{:,.4f}».format(np.percentile(V, 30))
P30 = «{:.0%}».format(0.30)
Per35 = «{:,.4f}».format(np.percentile(V, 35))
P35 = «{:.0%}».format(0.35)
Per40 = «{:,.4f}».format(np.percentile(V, 40))
P40 = «{:.0%}».format(0.40)
Per45 = «{:,.4f}».format(np.percentile(V, 45))
P45 = «{:.0%}».format(0.45)
Per50 = «{:,.4f}».format(np.percentile(V, 50))
P50 = «{:.0%}».format(0.50)
Per55 = «{:,.4f}».format(np.percentile(V, 55))
P55 = «{:.0%}».format(0.55)
Per60 = «{:,.4f}».format(np.percentile(V, 60))
P60 = «{:.0%}».format(0.60)
Per65 = «{:,.4f}».format(np.percentile(V, 65))
P65 = «{:.0%}».format(0.65)
Per70 = «{:,.4f}».format(np.percentile(V, 70))
P70 = «{:.0%}».format(0.70)
Per75 = «{:,.4f}».format(np.percentile(V, 75))
P75 = «{:.0%}».format(0.75)
Per80 = «{:,.4f}».format(np.percentile(V, 80))
P80 = «{:.0%}».format(0.80)
Per85 = «{:,.4f}».format(np.percentile(V, 85))
P85 = «{:.0%}».format(0.85)
Per90 = «{:,.4f}».format(np.percentile(V, 90))
P90 = «{:.0%}».format(0.90)
Per95 = «{:,.4f}».format(np.percentile(V, 95))
P95 = «{:.0%}».format(0.95)
Minimum = «{:,.4f}».format(np.min(V))
Maximum = «{:,.4f}».format(np.max(V))
Range = «{:,.4f}».format((np.max(V))-(np.min(V)))
Count = len(V)
Mean = «{:,.4f}».format(np.mean(V))
Median = «{:,.4f}».format(np.median(V))
Mode = «{:,.4f}».format(stats.mode(V)[0][0])
Variance = np.var(V)
Std_Dev = «{:,.4f}».format(np.std(V))
Std_Err = «{:,.4f}».format(np.std(V)/np.sqrt(Count))
Skewness = round(stats.skew(V),9)
Kurtosis = round((stats.kurtosis(V)+3),9)
Exc_Kurt = round((stats.kurtosis(V)),9)
Left_X = Per5
Left_P = P5
Right_X = Per95
Right_P = P95
Diff_X = «{:,.4f}».format((np.percentile(V, 95) — np.percentile(V, 5)))
Diff_P = «{:.0%}».format(0.90)
g = {«Information»: e, «Result»: f}
st = pd.DataFrame(data=g)
a = [«Minimum»,»Maximum»,»Range»,»Count»,»Mean»,»Median»,»Mode»,»Variance»,»Std Dev»,»Std Err»,»Skewness»,»Kurtosis»\
,»Exc_Kurt»,»Left X»,»Left P»,»Right X»,»Right P»,»Diff X»,»Diff P»]
b = [Minimum, Maximum, Range, Count, Mean, Median, Mode, Variance, Std_Dev, Std_Err, Skewness, Kurtosis\
,Exc_Kurt, Left_X, Left_P, Right_X ,Right_P, Diff_X, Diff_P]
c = [P5,P10,P15,P20,P25,P30,P35,P40,P45,P50,P55,P60,P65,P70,P75,P80,P85,P90,P95]
d = [Per5, Per10, Per15, Per20, Per25, Per30, Per35, Per40, Per45, Per50, Per55, Per60, Per65,\
Per70, Per75, Per80, Per85, Per90, Per95]
d = {«Statistics»: a, «Statistics Result»: b, «Percentile»: c, «Percentile Result»: d}
st1 = pd.DataFrame(data=d)
from datetime import date
today = date.today()
now = datetime.now()
import calendar
curr_date = date.today()
print(«\033[1m Simulation Summary Information»)
print(«\033[0m ================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))
st

В-седьмых, проанализируем статистику будущей стоимости собственного капиталаprint(«\033[1m Summary Statistics for Future Value of Equity»)
print(«\033[0m ======================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))
st1

В-восьмых, давайте визуализируем функцию плотности вероятности (pdf) для будущей стоимости собственного капиталаplt.figure(figsize = (4,4))
sns.set(font_scale = 1.2)
ax = sns.histplot(data=V,bins=500,color=»red»)
ax.set_xlabel(«Values»,fontsize=14)
ax.set_ylabel(«Frequency»,fontsize=14)
ax.set_xlim(np.percentile(V, 2.5),np.percentile(V, 97.5))
print(«\033[1m Probability Density Function for the Future Value of Equity (Sim#1)»)
print(«\033[0m ===================================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))

Итак, давайте визуализируем функцию кумулятивного распределения (cdf) для будущей стоимости собственного капиталаplt.figure(figsize = (4,4))
kwargs = {«cumulative»: True}
sns.axes_style(«whitegrid»)
sns.set(font_scale = 1.2)
ax = sns.distplot(V, hist_kws=None, kde_kws=kwargs,color=»red»)
ax.set_xlabel(«Values»,fontsize=14)
ax.set_ylabel(«Frequency»,fontsize=14)
ax.set_xlim(np.percentile(V, 2.5),np.percentile(V, 97.5))
print(«\033[1m Cumulative Distribution Function for the Future Value of Equity (Sim#1)»)
print(«\033[0m =======================================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))

Теперь найдем оценку будущей стоимости собственного капиталаV = df[‘ST’]
print(«The future value of equity is: {:.2f}».format(np.mean(V)))

Будущая стоимость собственного капитала составляет: 211,32

Давайте оценим точность нашей оценкиprint(«The stochastic error of the future value of equity is: {:.2f}».format(np.std(V)/np.sqrt(Paths)))

Стохастическая ошибка будущей стоимости эквити составляет: 1,44

Стандартные отклонения нормы, соответствующие 95% доверительному интервалу, составляют αMIN = −1,96 и αMAX = 1,96. Другими словами, у нас есть 2,5% в каждом хвосте.print(«The lower bound for the future value of equity with 95% confidence level is {:.4}».format(np.mean(V) — 1.96*np.std(V)/np.sqrt(Paths)))

Нижняя граница будущей стоимости акций с доверительной вероятностью 95% составляет 208,4932print(«The upper bound for the future value of equity with 95% confidence level is {:.4}».format(np.mean(V) + 1.96*np.std(V)/np.sqrt(Paths)))

Верхняя граница будущей стоимости акций с доверительной вероятностью 95% составляет 214,1442

8 Заключение

Мы видим, что наша стохастическая мировая оценка будущей стоимости капитала ($211,3187) очень близка к нашей детерминированной мировой оценке будущей стоимости капитала ($211,7000). Это означает, что мы строим хорошую модель для прогнозирования будущей стоимости капитала, используя прогнозное моделирование по методу Монте-Карло.

Построение вероятностной оценки риска с использованием моделирования по методу Монте-Карло с помощью Python MCerp

В этой статье объясняется, как создать стохастическую оценку риска с помощью моделирования по методу Монте-Карло с использованием библиотеки Python MCerp

Источник изображения: Адаптировано из PlugPng

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

https://cdn.embedly.com/widgets/media.html?src=https%3A%2F%2Fwww.youtube.com%2Fembed%2FUISccjS1Gv4%3Ffeature%3Doembed&display_name=YouTube&url=https%3A%2F%2Fwww.youtube.com%2Fwatch%3Fv%3DUISccjS1Gv4&image=https%3A%2F%2Fi.ytimg.com%2Fvi%2FUISccjS1Gv4%2Fhqdefault.jpg&key=a19fcc184b9711e1b4764040d3dc5c07&type=text%2Fhtml&schema=youtube

Если вы столкнулись с этой проблемой, я бы посоветовал прочитать оригинальную статью 2020 года ниже, поскольку мы воссоздадим ту же вероятностную оценку праздничного бюджета, но на этот раз с использованием Python вместо надстройки Excel.

Если вы хотите следовать за нами, Jupyter Notebook доступен в этом репозитории Github.

GitHub — ZhijingEu/MonteCarloSim-Via-PythonMCerp-Example: Это простой пример того, как запустить…

Это простой пример того, как запустить моделирование по методу Монте-Карло для разработки анализа рисков бюджетных затрат с помощью Python…

github.com

Очертание

1.Краткий обзор библиотек моделирования по методу Монте-Карло в Python
2.Тематическое исследование Фоновый остров фантазий (Redux)
3. Настройка входных и выходных переменных
4.Запуск симуляции и агрегирование результатов
5.Моделирование зависимостей между неопределенными переменными — корреляция
6.Моделирование зависимостей — мультипликативные «метапеременные»
7.Анализ ключевых факторов, влияющих на результаты
8.Заключение

1.Краткий обзор библиотек моделирования по методу Монте-Карло в Python

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

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

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

Код, на который ссылается ниже, в значительной степени заимствован из этой статьи другого автора Medium Хайко Оннена (Heiko Onnen), у которого также есть ряд других статей, связанных с более продвинутыми темами по стохастическому моделированию.

2. Предыстория кейса: Остров фантазий (Redux)

В качестве краткого резюме — пример, который мы будем анализировать, включает в себя упражнение по планированию бюджета для отпуска в вымышленной зимней островной стране, чтобы ответить на вопрос «Сколько будет стоить мой отпуск на острове фантазий?».

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

Вот! Остров фантазий — место отдыха вашей мечты (хотя, возможно, вы захотите пропустить рафтинг)

3. Настройка входных и выходных переменных

Входные данные для вероятностных оценок обычно состоят из:

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

Для начала нам сначала нужно импортировать пару библиотек%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from mcerp import correlate, correlation_matrix, plotcorr
from mcerp import Beta, Uniform, PERT, Binomial
from mcerp import uv, stats

from scipy import stats as stats
from scipy.stats import rv_continuous, rv_histogram

np.set_printoptions(precision=3, suppress=True)
pd.options.display.float_format = ‘{:.3f}’.format

MCerp предоставляет ряд различных распределений вероятностей, которые построены с помощью scipy.stats, но имеют несколько полезных функций. В качестве примера, вот как настроить дистрибутив Beta Pert и что MCerp предоставляет в качестве выходных данныхnpts=10000 #No of simulation runs

x=PERT(1,4,10)
x.describe()
rx=x._mcpts # shows an array of all the individual simulation results
rx

Мы также создадим вспомогательную функцию для построения графика функции распределения вероятностейdef plot_hist(data, title):
hist = np.histogram(data, bins=100)
histdist = rv_histogram(hist)

X = np.linspace(data.min(), data.max(), 100)
plt.title(title)
plt.hist(data, density=True, bins=100)
plt.plot(X, histdist.pdf(X))
plt.show()

rx=x._mcpts
plot_hist(rx,»Example PERT Distribution»)

Теперь, когда все это настроено, давайте определим все входные переменные, как показано в примере

#InputVariables

PlaneFare= PERT(5000,7000,12000,PertLambda) #BetaPert
AccommodationCosts = PERT(2000,2500,4000,PertLambda) #BetaPert
Meals=PERT(500,600,900,PertLambda) #BetaPert
MiscShoppingExpenses=1000 #Has No Range as It is UnRisked
HolidayActivityTours=Uniform(1500,4000) #Uniform
ClothingTravelGear=PERT(300,400,500,PertLambda) #BetaPert

MedicalEmergencyRiskLikehood=Binomial(1, 0.1)
MedicalEmergencyRiskImpact=Uniform(1000,5000)

TheftRiskLikelihood=Binomial(1, 0.1)
TheftRiskImpact=PERT(500,750,2000,PertLambda) #BetaPert#OutputVariables
RiskedRanges = PlaneFare+AccommodationCosts+Meals+MiscShoppingExpenses+HolidayActivityTours+ClothingTravelGear
DiscreteRisks=MedicalEmergencyRiskLikehood*MedicalEmergencyRiskImpact+TheftRiskLikelihood*TheftRiskImpact
TotalHolidayCosts=RiskedRanges+DiscreteRisks

4.Запуск симуляции и агрегирование результатов

Нам по-прежнему нужно собрать все отдельные наборы результатов в связный табличный формат, что и делается в следующих нескольких строках кода, где я использовал префикс «r» для обозначения массива, в котором хранятся все значения длясоответствующей входной или выходной переменной# collect the results in array variables

rPlaneFare = PlaneFare._mcpts
rAccommodationCosts = AccommodationCosts._mcpts
rMeals = Meals._mcpts
rHolidayActivityTours = HolidayActivityTours._mcpts
rClothingTravelGear = ClothingTravelGear._mcpts
rMedicalEmergencyRiskLikehood = MedicalEmergencyRiskLikehood._mcpts
rMedicalEmergencyRiskImpact=MedicalEmergencyRiskImpact._mcpts
rTheftRiskLikelihood=TheftRiskLikelihood._mcpts
rTheftRiskImpact=TheftRiskImpact._mcpts
rRiskedRanges=RiskedRanges._mcpts
rDiscreteRisks=DiscreteRisks._mcpts
rTotalHolidayCosts=TotalHolidayCosts._mcpts

# combine the arrays in a 2-dimensional array
rand1 = np.vstack((rPlaneFare,rAccommodationCosts,rMeals,rHolidayActivityTours,rClothingTravelGear,
rMedicalEmergencyRiskLikehood,rMedicalEmergencyRiskImpact,
rTheftRiskLikelihood,rTheftRiskImpact,rRiskedRanges,rDiscreteRisks,rTotalHolidayCosts))

VarNames=[«Input_PlaneFare», «Input_AccomCosts», «Input_Meals», «Input_HolidayActivityTours», «Input_ClothingTravelGear»,
«Input_MedEmergencyChance»,»Input_MedEmergencyCost»,»Input_TheftRiskChance»,»Input_TheftRiskCost»,
«Output_RiskedRanges»,»Output_DiscreteRisks»,»Output_TotalHolidayCosts»]

# copy the array to a dataframe for a more transparent layout
df1 = pd.DataFrame(data=rand1).T
df1.columns=VarNames
df1

Используя метод describe в MCerp, мы можем получить некоторую ключевую статистику

Мы также вычислим 10-й, 25-й, медианный (50-й), 75-й, 90-й процентили из 10 000 запусков моделирования с использованием метода квантиля пандыdf1.quantile([.1, .25, .5, .75, .9], axis = 0)

Однако, поскольку их немного сложно визуализировать, давайте создадим несколько графиков функции распределения вероятностейVarSimResults=[rPlaneFare,rAccommodationCosts,rMeals,rHolidayActivityTours,rClothingTravelGear,
rMedicalEmergencyRiskLikehood,rMedicalEmergencyRiskImpact,rTheftRiskLikelihood,
rTheftRiskImpact,rRiskedRanges,rDiscreteRisks,rTotalHolidayCosts]

for n in range(0,len(VarSimResults)):
plot_hist(VarSimResults[n],VarNames[n])

5.Моделирование зависимостей между неопределенными переменными — корреляция

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

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

Например, построение графика (некоррелированного) рискованного диапазона входных данных# get the correlation matrix BEFORE applying correlation
BeforeCorrApplied = correlation_matrix([AccommodationCosts, Meals, HolidayActivityTours])
BeforeCorrApplied

# plot correlations BEFORE applying correlation
corrplot0 = plotcorr([AccommodationCosts, Meals, HolidayActivityTours], labels=[«Input_Accom»,
«Input_Meals», «Input_HolidayActivityTours»])

Значения, близкие к 0,00, указывают на то, что корреляция между выбранными входными переменными незначительна или отсутствует

Если мы уже определили, что существует определенное корреляционное поведение, которое мы ожидаем между переменными, у MCerp есть функция, которая может применить эту корреляцию с помощью математического метода, называемого Иман-Коновером, который приблизительно индуцирует заданный ранговый порядок (т.е. корреляцию Спирмена) между переменными-компонентами, не изменяя одномерного распределения компонентов

(Вот статья из еще более продвинутой библиотеки Python под названием aggregate, которую мы не применяем в примере, но все же дает довольно подробное объяснение того, как она работает)

aggregate.readthedocs.ioapplied_correlation=np.array([[1.00, 0.85, 0.85],
[0.85, 1.00, 0.85],
[0.85, 0.85, 1.0]])

# impose the targeted correlation matrix on the 3 input variables
correlate([AccommodationCosts, Meals, HolidayActivityTours], applied_correlation)

# plot the new correlation matrix of the input variables
corrplot2 = plotcorr([AccommodationCosts, Meals, HolidayActivityTours],
labels=[«Input_Accom»,»Input_Meals», «Input_HolidayActivityTours»])

Следующим шагом является подтверждение того, что корреляции были индуцированы до значения, достаточно близкого к ожидаемым целевым корреляциям

0.83 и 0.846 достаточно близко к 0.85, я полагаю

Затем нам нужно «сбросить» выходные переменные, чтобы увидеть эффект#OutputVariables — Rerun with Corr Variables
RiskedRanges = PlaneFare+AccommodationCosts+Meals+MiscShoppingExpenses+HolidayActivityTours+ClothingTravelGear
DiscreteRisks=MedicalEmergencyRiskLikehood*MedicalEmergencyRiskImpact+TheftRiskLikelihood*TheftRiskImpact
TotalHolidayCosts=RiskedRanges+DiscreteRisks

и повторно запустите симуляцию# run the simulation model,
# and collect the results in array variables

rPlaneFare = PlaneFare._mcpts
rAccommodationCosts = AccommodationCosts._mcpts
rMeals = Meals._mcpts
rHolidayActivityTours = HolidayActivityTours._mcpts
rClothingTravelGear = ClothingTravelGear._mcpts
rMedicalEmergencyRiskLikehood = MedicalEmergencyRiskLikehood._mcpts
rMedicalEmergencyRiskImpact=MedicalEmergencyRiskImpact._mcpts
rTheftRiskLikelihood=TheftRiskLikelihood._mcpts
rTheftRiskImpact=TheftRiskImpact._mcpts
rRiskedRanges=RiskedRanges._mcpts
rDiscreteRisks=DiscreteRisks._mcpts
rTotalHolidayCosts=TotalHolidayCosts._mcpts

# combine the arrays in a 2-dimensional array
rand1 = np.vstack((rPlaneFare,rAccommodationCosts,rMeals,rHolidayActivityTours,rClothingTravelGear,
rMedicalEmergencyRiskLikehood,rMedicalEmergencyRiskImpact,
rTheftRiskLikelihood,rTheftRiskImpact,rRiskedRanges,rDiscreteRisks,rTotalHolidayCosts))

VarNames=[«Input_PlaneFare», «Input_AccomCosts», «Input_Meals», «Input_HolidayActivityTours», «Input_ClothingTravelGear»,
«Input_MedEmergencyChance»,»Input_MedEmergencyCost»,»Input_TheftRiskChance»,»Input_TheftRiskCost»,
«Output_RiskedRanges»,»Output_DiscreteRisks»,»Output_TotalHolidayCosts»]

# copy the array to a dataframe for a more transparent layout
df_postCorr = pd.DataFrame(data=rand1).T
df_postCorr.columns=VarNames
df_postCorr

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

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

По причинам, которые станут ясны в следующем разделе, давайте также конвертируем бюджет, деноминированный в Fantasy Bucks, в доллары США, используя фиксированное соотношение 1:4

Среднее значение в долларах США — это просто среднее значение фэнтезийных баксов/4, а дисперсия в долларах США — это дисперсия в фэнтези-баксах/4²

6.Моделирование зависимостей — мультипликативные «метапеременные»

Все шаги, описанные выше, отображают результаты в виде фэнтезийных баксов, но в этом вымышленном примере нам нужно планировать наш бюджет в долларах США.

Чтобы обеспечить большую вариативность результата, мы также можем смоделировать неопределенность FOREX, которая «усиливает» неопределенность общего результата (поскольку другие позиции, которые сами по себе имеют неопределенность, умножаются на этот переменный курс FOREX)FOREXVariability=PERT(3.75,4.00,4.50) #BetaPert

#OutputVariables — Rerun with Corr Variables
RiskedRanges = PlaneFare+AccommodationCosts+Meals+MiscShoppingExpenses+HolidayActivityTours+ClothingTravelGear
DiscreteRisks=MedicalEmergencyRiskLikehood*MedicalEmergencyRiskImpact+TheftRiskLikelihood*TheftRiskImpact
TotalHolidayCosts=RiskedRanges+DiscreteRisks
TotalHolidayCostsInUSD=TotalHolidayCosts/FOREXVariability

Предполагая, что предыдущие результаты были с фиксированной ставкой 4 Fantasy Bucks = 1 USD, и сравнивая результаты, изменчивость FOREX создает более широкий диапазон от P10 до P90

7.Анализ ключевых факторов, влияющих на результаты

Как и в предыдущем примере XLRisk, после того, как у вас есть выходные данные, ключевым моментом является определение того, какие входные данные риска являются основными факторами конечного выходного значения, которое нас интересует#Viewing the results by ranking the correlations between each random input variable
#(i.e uncertain ranges or risk events likelihoods or impact ranges) and the outcome

# check the new correlation matrix for the 3 input variables
OutputVsInputCorr = correlation_matrix([TotalHolidayCostsInUSD,PlaneFare,AccommodationCosts,
Meals,HolidayActivityTours,ClothingTravelGear,
MedicalEmergencyRiskLikehood, MedicalEmergencyRiskImpact,
TheftRiskLikelihood,TheftRiskImpact,
FOREXVariability])

df_Correlation=pd.DataFrame(data=OutputVsInputCorr)
df_Correlation.columns=[«Output_TotalHolidayCosts»,»Input_PlaneFare», «Input_AccomCosts», «Input_Meals»,
«Input_HolidayActivityTours», «Input_ClothingTravelGear»,
«Input_MedEmergencyChance»,»Input_MedEmergencyCost»,»Input_TheftRiskChance»,»Input_TheftRiskCost»,
«FOREX Variability»]

df_Correlation.T[0][1:].sort_values(ascending=False)

df_Correlation.T[0][1:].sort_values().plot.barh()

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

8.Заключение

1. Моделирование по методу Монте-Карло работает путем многократной выборки из набора заранее определенных случайных величин в модели и агрегирования общих результатов для понимания общего диапазона результатов

2.Исходные данные для вероятностных оценок обычно состоят из:

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

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

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

Надеемся, что эта статья послужила для вас хорошим введением в использование моделирования по методу Монте-Карло для построения вероятностных оценок риска.

9. Изучение других генераторов случайных чисел и методов выборки

Существует несколько различных методов оптимизации «подбрасывания монеты», которое происходит для каждой итерации в симуляции по методу Монте-Карло. Одна из областей связана с тем, как мы выполняем случайную выборку для различных распределений

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

Наиболее распространенным методом, используемым в коммерческом программном обеспечении для моделирования (Palisade @Risk, Vose Software Model Risk, Oracle Crystal Ball), является так называемая латинская гиперкубическая выборка, при которой каждая переменная разбивается на равномерно распределенные области, из которых выбираются случайные выборки для получения контролируемой случайной выборки.

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

К сожалению, MCerp по умолчанию использует только Latin Hypercube Sampling и не позволяет переключаться на другие методы сэмплирования. Тем не менее, другие библиотеки Python, такие как Monaco, имеют эту опцию, и если вы хотите, вы даже можете сделать это с помощью базового движка Scipy Quasi Monte Carlo, который позволяет вам переключаться между Latin Hypercube, Sobol или Halton Sampling.

Помимо методологии случайной выборки, еще одной областью, с которой мы могли бы поиграть, является функция генератора псевдослучайных чисел (PRNG). В большинстве коммерческих инструментов моделирования ГПСЧ, как правило, представляет собой нечто, называемое мерсеннским твистером, которому отдают предпочтение из-за его вычислительной эффективности

К сожалению, я не смог понять, как это изменить в MCerp. Похоже, что MCerp и все другие библиотеки, на которые ссылаются (Monaco и PyMCSL), построены поверх SciPy, который, согласно документации, основан на Numpy, который использует что-то еще, называемое Permuted Congruential Generator (64-bit, PCG64) в качестве PRNG по умолчанию

Все читатели, которым удалось дойти до этого места в статье и могут быть более осведомленными — пожалуйста, оставьте комментарий, если вы знаете, как переключить PRNG по умолчанию в Numpy (и стоит ли это делать…)

10.За пределами корреляций с связками

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

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

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

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

Математика действительно задействована, но интуиция, которую мы можем вынести, заключается в том, что метод Имана Коновера аналогичен применению гауссовской связки, где «индуцированная связь» между переменными одинакова во всем диапазоне возможных значений.

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

Чтобы вы не подумали, что это неясная и эзотерическая тема, обратите внимание на следующую поучительную историю…

Когда в 2000 году была впервые опубликована функция Гаусса Дэвида Ли, инвесторы использовали ее как быстрый — и фатально ошибочный — способ оценить сотни миллиардов долларов CDO, заполненных ипотечными кредитами. К сожалению, в 2008 году, когда рынки начали вести себя так, что модель не могла с этим справиться, это привело к потере триллионов долларов и поставило под серьезную угрозу выживание глобальной банковской системы.

11.Рекомендации по проведению вероятностного анализа

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

11.1 Обрамление модели и допущения

Источник Изображения: Alamy
  • Почему выбор моделируется как неопределенность?
  • Например, если диапазоны неопределенности стоимости отражают спектр от пятизвездочного отеля до молодежного хостела, то, возможно, лучшим подходом может быть предложение более простых и менее детализированных сценариев и их влияния на стоимость, т.е. Роскошный маршрут отдыха против маршрута туриста с завязкой обуви
  • Таким образом, более детальный вероятностный анализ бюджета следует проводить только ПОСЛЕ того, как первоначальное решение будет сужено

11.2 Последовательность в сборе и определении рисков

Источник Изображения: Blogspot
  • Насколько мы уверены в том, что диапазоны неопределенности и значения влияния рисковых событий не удваиваются и применяются последовательно?
  • В приведенном выше примере есть диапазон одежды/снаряжения для путешествий, и мы обозначили его в фэнтезийных баксах, но если мы планируем купить его перед поездкой в интернет-магазине в долларах США, неверно применять неопределенность FOREX к этой позиции.
  • Это те вещи, на которые следует обращать внимание, особенно когда модель включает неопределенности «второго порядка», которые имеют мультипликативный эффект
  • Например, в анализе рисков по расписанию существует популярная методика, называемая методом драйверов риска, которая позволяет аналитикам избежать необходимости беспокоиться о взаимозависимостях/корреляциях, но предполагает, что мы можем оценить % влияния hi-mid-low взаимоисключающим образом (что на практике практически невозможно оценить)
  • Кроме того, если модель включает в себя «условные риски» (т.е. события риска для детей, которые имеют шанс произойти только в том случае, если произойдет родительское событие риска), вам нужно будет откалибровать оценки вероятности (в качестве тривиального примера, вероятность того, что случайно выбранный человек заболеет раком легких, была бы другой, если бы мы знали, что случайный человек является заядлым курильщиком)

11.3 Преобразование полученных результатов в ДЕЙСТВЕННЫЕ выводы

  • Компетентный аналитик должен не только объяснить основы анализа, модельные предположения и выводы, но и сделать полученные выводы действенными.
  • Например, недостаточно подчеркнуть, что изменчивость стоимости авиабилетов оказывает наибольшее влияние на общую неопределенность бюджета, но также проверить базовые допущения модели и понять факторы, влияющие на стоимость отдыхающих:
  • >>Если отпуск ограничен в бюджете, то, возможно, отдыхающим нужно согласиться на рейсы в «нечетное время», чтобы снизить риск затрат и подчеркнуть, что чем дольше они откладывают бронирование рейса, тем дороже он становится (т.е. обеспечение того, чтобы решения были ОГРАНИЧЕНЫ ВРЕМЕНЕМ)
  • >>Еще более серьезный вопрос заключается в том, чтобы понять гибкость дат — т.е. можно ли перенести время отпуска на более позднюю дату, когда он не пик и Остров Фантазий не наводнен туристами, чтобы насладиться более низкими ценами (но с учетом того, что некоторые достопримечательности могут быть закрыты)

12. В заключение… (На этот раз по-настоящему)

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

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

Моделирование по методу Монте-Карло и треугольное распределение в Python

Треугольное распределение описывает ситуацию, в которой известны минимальные, максимальные и наиболее вероятные значения. Например, можно описать объем продаж продукта, когда прошлые данные об объеме продаж показывают минимальный, максимальный и наиболее вероятный объем продаж этого продукта.

Фото: Внутренняя стоимость

Мотивация

Чтобы спрогнозировать объем продаж в следующем году, я решил использовать моделирование по методу Монте-Карло и треугольное распределение. Моделирование по методу Монте-Карло — это удивительный способ позаботиться о неопределенности в наших моделях.

Допустим, мы хотим спрогнозировать объем продаж нашей компании в следующем году, скажем, на 2023 год. Предполагая, что прошлые данные об объеме продаж показывают, что минимальный объем продаж составил 25 000 долларов США, максимальный объем продаж — 45 000 долларов США, а объем продаж в режиме — 40 000 долларов США, как мы можем использовать моделирование по методу Монте-Карло для прогнозирования объема продаж на 2023 год на основе этих данных?Minimum = 25000
MostLikely = 40000
Maximum = 45000

Итак, первое, что нам нужно сделать, это предположить распределение по объему продаж. Так как здесь мы имеем минимум, наиболее вероятный (т.е. моду) и максимум — значит, это треугольное распределение.

Моделирование по методу Монте-Карло

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

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

Блестящий аналитик будет использовать все доступные инструменты, чтобы получить тот же ответ самым простым и практичным способом. И во всех случаях, при правильном моделировании, моделирование по методу Монте-Карло дает те же ответы, что и более математически элегантные методы.

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

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

Упрощенно можно представить себе подход моделирования по методу Монте-Карло как многократное извлечение мячей для гольфа из большой корзины с заменой. Размер и форма корзины зависят от предположения о распределении входных данных (например, нормальное распределение со средним значением 100 и стандартным отклонением 10, по сравнению с равномерным распределением или треугольным распределением), когда некоторые корзины глубже или более симметричны, чем другие, что позволяет вытаскивать определенные шары чаще, чем другие.

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

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

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

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

Треугольное распределение

Треугольное распределение описывает ситуацию, в которой известны минимальные, максимальные и наиболее вероятные значения. Например, можно описать объем продаж продукта, когда прошлые данные об объеме продаж показывают минимальный, максимальный и наиболее вероятный объем продаж этого продукта.

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

Условия

Три условия, лежащие в основе треугольного распределения:

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

Математические конструкции для треугольного распределения следующие:

Минимальное значение (Min), наиболее вероятное значение (Likely) и максимальное значение (Max) являются параметрами распределения.

Требования к входным данным:

Min ≤ Most Likely ≤ Max и может принимать любое значение.

Однако Min < Max и может принимать любое значение.

Предварительная обработка данных

Прежде всего, давайте импортируем все соответствующие библиотеки:import numpy as np
import pandas as pd
from scipy.stats import norm
from numpy.random import randn
from numpy import random as rn
import scipy.stats as si
from matplotlib import pyplot as plt
import seaborn as sns
from matplotlib import pyplot as plt
from IPython.display import Image
%matplotlib inline

Во-вторых, давайте создадим фрейм данных из 50 000 симуляций (т.е. испытаний) случайных чисел, которые равномерно распределены между 0 и 1.Simulations = 50000
f = np.arange(1,Simulations+1)
df = pd.DataFrame(data=f)
df[‘RAND’] = rn.rand(Simulations,1)
df.drop([0], axis=1, inplace=True)
df

И тогда для треугольного распределения есть формула

У нас есть более низкий диапазон, более высокий диапазон, у нас есть наш общий диапазон и наша кумулятивная вероятность, которая на самом деле является нашими случайными числами (т.е. df["RAND"]). Итак, у нас есть генерация случайных чисел с треугольным распределением.

Давайте просто отформатируем и их. Наш нижний диапазон — это, скорее всего, наш минус наш минимум. Наш верхний диапазон — это наш максимум минус наш наиболее вероятный диапазон, а наш общий диапазон — это наш максимум минус наш минимум.MostLikely = Mode
LowerRange = Mode — Minimum
HigherRange = Maximum — Mode
TotalRange = Maximum — Minimum

Приведенная выше формула немного сложна для понимания, но давайте просто введем ее. Таким образом, нам нужно проверить, находится ли кумулятивная вероятность (т.е. df["RAND"][i]«][i] ) ниже нижнего диапазона над общим диапазоном, а затем наш объем продаж (т.е. df["SalesVolume"][i] ) равен нашему минимуму плюс квадратный корень из нашей кумулятивной вероятности, умноженный на наш нижний диапазон, умноженный на наш общий диапазон, и если нет, то нет (т.е. else), то наш объем продаж равен нашему максимальному минус квадратный корень из 1 минус наша кумулятивная вероятность, умноженная на больший диапазон, умноженный на наш общий диапазон.

И самый простой способ проверить, правильно ли мы написали нашу формулу, заключается в том, что если наша кумулятивная вероятность равна 0, нам нужно получить 25 000, что является минимумом, а если она равна 1, нам нужно получить 45 000, что является максимумом. Итак, наша формула верна.df[‘SalesVolume’] = 0
for i in range(0,len(df[‘SalesVolume’])):
if df[‘RAND’][i] < (LowerRange)/(TotalRange):
df[‘SalesVolume’][i] = Minimum + np.sqrt(df[‘RAND’][i]*(LowerRange)*(TotalRange))
else:
df[‘SalesVolume’][i] = Maximum — np.sqrt((1-df[‘RAND’][i])*(HigherRange)*(TotalRange))
df

Обозначим наш целевой столбец буквой V.V = df[‘SalesVolume’]
V

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

Входными данными функции является вероятность, представляющая область под стандартной нормальной кумулятивной функцией распределения (cdf) в диапазоне от минус бесконечности до плюс бесконечности. Выходными данными функции является число стандартных отклонений, соответствующих данной конкретной области при стандартной нормальной кумулятивной функции распределения.def NORMSINV(x):
x = si.norm.ppf(x)
return (x)

Например, обратное стандартному нормальному кумулятивному распределению, с вероятностью 0,9500 равно 1,6448536269514722NORMSINV(0.95)

1.6448536269514722

Анализ данных

Проанализируем информацию моделированияG = 857550884
np.random.seed(G)
from datetime import datetime
start_time = datetime.now()
from datetime import datetime
# datetime object containing current date and time
now = datetime.now()
dt_string = now.strftime(«%d/%m/%Y %H:%M:%S»)

time_elapsed1 = datetime.now() — start_time
from scipy import stats
Workbook_Name = «Triangular Distribution with MC Simulations.ipynb»
Number_of_Simulations = «{:,.0f}».format(Simulations)
Number_of_Iterations = «{:,.0f}».format(1)
Number_of_Inputs = «{:,.0f}».format(3)
Number_of_Outputs = 38
Sampling_Type = «Latin Hypercube»
Simulation_Start_Time = dt_string
Simulation_Duration = «{}».format(time_elapsed1)
Random_N_Generator = «Mersenne Twister»
Random_Seed = G
e = [«Workbook Name»,»Number of Simulations»,»Number of Iterations»,»Number of Inputs»,»Number of Outputs»,»Sampling Type»,\
«Simulation Start Time»,»Simulation Duration»,»Random # Generator»,»Random Seed»]
f = [Workbook_Name, Number_of_Simulations, Number_of_Iterations, Number_of_Inputs, Number_of_Outputs, Sampling_Type,\
Simulation_Start_Time, Simulation_Duration, Random_N_Generator, Random_Seed]
Per5 = «{:,.1f}».format(np.percentile(V, 5))
P5 = «{:.0%}».format(0.05)
Per10 = «{:,.1f}».format(np.percentile(V, 10))
P10 = «{:.0%}».format(0.10)
Per15 = «{:,.1f}».format(np.percentile(V, 15))
P15 = «{:.0%}».format(0.15)
Per20 = «{:,.1f}».format(np.percentile(V, 20))
P20 = «{:.0%}».format(0.20)
Per25 = «{:,.1f}».format(np.percentile(V, 25))
P25 = «{:.0%}».format(0.25)
Per30 = «{:,.1f}».format(np.percentile(V, 30))
P30 = «{:.0%}».format(0.30)
Per35 = «{:,.1f}».format(np.percentile(V, 35))
P35 = «{:.0%}».format(0.35)
Per40 = «{:,.1f}».format(np.percentile(V, 40))
P40 = «{:.0%}».format(0.40)
Per45 = «{:,.1f}».format(np.percentile(V, 45))
P45 = «{:.0%}».format(0.45)
Per50 = «{:,.1f}».format(np.percentile(V, 50))
P50 = «{:.0%}».format(0.50)
Per55 = «{:,.1f}».format(np.percentile(V, 55))
P55 = «{:.0%}».format(0.55)
Per60 = «{:,.1f}».format(np.percentile(V, 60))
P60 = «{:.0%}».format(0.60)
Per65 = «{:,.1f}».format(np.percentile(V, 65))
P65 = «{:.0%}».format(0.65)
Per70 = «{:,.1f}».format(np.percentile(V, 70))
P70 = «{:.0%}».format(0.70)
Per75 = «{:,.1f}».format(np.percentile(V, 75))
P75 = «{:.0%}».format(0.75)
Per80 = «{:,.1f}».format(np.percentile(V, 80))
P80 = «{:.0%}».format(0.80)
Per85 = «{:,.1f}».format(np.percentile(V, 85))
P85 = «{:.0%}».format(0.85)
Per90 = «{:,.1f}».format(np.percentile(V, 90))
P90 = «{:.0%}».format(0.90)
Per95 = «{:,.1f}».format(np.percentile(V, 95))
P95 = «{:.0%}».format(0.95)
Minimum = «{:,.1f}».format(np.min(V))
Maximum = «{:,.1f}».format(np.max(V))
Mean = «{:,.1f}».format(np.mean(V))
Std_Dev = «{:,.1f}».format(np.std(V))
Variance = np.var(V)
Std_Error = «{:,.1f}».format(np.std(V)/np.sqrt(Simulations))
Skewness = round(stats.skew(V),9)
Kurtosis = round((stats.kurtosis(V)+3),9)
Median = «{:,.1f}».format(np.median(V))
Mode = «{:,.1f}».format(stats.mode(V)[0][0])
Left_X = Per5
Left_P = P5
Right_X = Per95
Right_P = P95
Diff_X = «{:,.1f}».format((np.percentile(V, 95) — np.percentile(V, 5)))
Diff_P = «{:.0%}».format(0.90)
Confidence_Interval = P95
Lower_Bound = «{:,.1f}».format((np.mean(V) — (np.std(V)/np.sqrt(Simulations))*NORMSINV(0.975)))
Upper_Bound = «{:,.1f}».format((np.mean(V) + (np.std(V)/np.sqrt(Simulations))*NORMSINV(0.975)))
g = {«Information»: e, «Result»: f}
st = pd.DataFrame(data=g)
a = [«Minimum»,»Maximum»,»Mean»,»Std Dev»,»Variance»,»Std Error», «Skewness»,»Kurtosis»,»Median»,»Mode»,\
«Left X»,»Left P»,»Right X»,»Right P»,»Diff X»,»Diff P»,»Confidence Interval»,»Lower Bound»,»Upper Bound»]
b = [Minimum, Maximum, Mean, Std_Dev, Variance, Std_Error, Skewness, Kurtosis, Median, Mode, Left_X, Left_P,\
Right_X, Right_P, Diff_X, Diff_P, Confidence_Interval, Lower_Bound, Upper_Bound]
c = [P5,P10,P15,P20,P25,P30,P35,P40,P45,P50,P55,P60,P65,P70,P75,P80,P85,P90,P95]
d = [Per5, Per10, Per15, Per20, Per25, Per30, Per35, Per40, Per45, Per50, Per55, Per60, Per65,\
Per70, Per75, Per80, Per85, Per90, Per95]
d = {«Statistics»: a, «Statistics Result»: b, «Percentile»: c, «Percentile Result»: d}
st1 = pd.DataFrame(data=d)
from datetime import date
today = date.today()
now = datetime.now()
import calendar
curr_date = date.today()
print(«\033[1m Simulation Summary Information»)
print(«\033[0m ================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))
st

Проанализируем статистику по прогнозируемому объему продажprint(«\033[1m Summary Statistics for Sales Volume»)
print(«\033[0m ======================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))
st1

Визуализация данных

Визуализируем функцию плотности вероятности (pdf) для прогнозируемого объема продажplt.figure(figsize = (4,4))
sns.set(font_scale = 1.2)
ax = sns.histplot(data=V,bins=10,color=»red»)
ax.set_xlabel(«Sales Volume»,fontsize=14)
ax.set_ylabel(«Frequency»,fontsize=14)
ax.set_xlim([20000,50000])
print(«\033[1m Probability Density Function for Sales Volume (Sim#1)»)
print(«\033[0m ======================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))

Визуализируем функцию кумулятивного распределения (cdf) для прогнозируемого объема продажplt.figure(figsize = (4,4))
kwargs = {«cumulative»: True}
sns.axes_style(«whitegrid»)
sns.set(font_scale = 1.2)
ax = sns.distplot(V, hist_kws=None, kde_kws=kwargs,color=»red»)
ax.set_xlabel(«Sales Volume»,fontsize=14)
ax.set_ylabel(«Frequency»,fontsize=14)
ax.set_xlim([20000,50000])
print(«\033[1m Cumulative Distribution Function for Sales Volume (Sim#1)»)
print(«\033[0m ==========================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))

Заключение

Для того, чтобы спрогнозировать объем продаж, я использовал моделирование по методу Монте-Карло и треугольное распределение. Моделирование методом Монте-Карло широко используется в ситуациях, когда отсутствует аналитическое решение и другие подходы (например, биномиальное дерево).

Моделирование по методу Монте-Карло использует моделирование факторов риска в соответствии с предполагаемым распределением вероятностей, а затем вычисляет объем продаж для каждого моделирования отдельно. Среднее значение объема продаж называется прогнозируемым значением объема продаж.

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

В прогнозах мы можем определить вероятность наступления, называемую доверительным интервалом. То есть, при двух значениях, каковы шансы, что результат окажется между этими двумя значениями?

В этих расчетах я использовал более 50 000 случайных симуляций. Существует 90% вероятность того, что окончательный результат (в данном случае объем продаж за 2023 год) составит от $28 916,9 до $42 759,0. С другой стороны, существует 5% вероятность того, что объем продаж будет ниже $28 916,9, и еще 5% вероятность того, что объем продаж будет выше $42 759,0.

То есть двусторонний доверительный интервал — это симметричный интервал, центрированный на среднем значении или медианном значении (т. е. значении 50-го процентиля). Таким образом, оба хвоста будут иметь одинаковую вероятность.

Среднее значение прогнозируемого объема продаж на 2023 год составляет $36 671,5. Определенно не значение в 40 000 долларов США в данных об объеме продаж за прошлые периоды и даже не среднее значение между минимальным значением данных об объеме продаж за прошлые периоды (т.е. 25 000 долларов США) и максимальным значением данных об объеме продаж за прошлые периоды (т.е. 45 000 долларов США), которое составляет 35 000 долларов США.

Моделирование по методу Монте-Карло и равномерное распределение в Python

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

Фото: Внутренняя стоимость

Мотивация

Чтобы спрогнозировать цену продажи в следующем году, я решил использовать моделирование по методу Монте-Карло и равномерное распределение. Моделирование по методу Монте-Карло — это удивительный способ позаботиться о неопределенности в наших моделях.

Допустим, мы хотим спрогнозировать цену продажи нашей компании в следующем году, скажем, на 2023 год. Предполагая, что прошлые данные о продажных ценах показывают, что минимальная цена продажи составила 8,00 долларов США, а максимальная цена продажи — 8,50 долларов США, как мы можем использовать моделирование по методу Монте-Карло для прогнозирования цены продажи на 2023 год на основе этих данных?Minimum = 8.00
Maximum = 8.50

Итак, первое, что нам нужно сделать, это предположить распределение цены продажи. Так как здесь мы имеем минимум и максимум — значит, это равномерное распределение.

Моделирование по методу Монте-Карло

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

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

Упрощенно можно представить себе подход моделирования по методу Монте-Карло как многократное извлечение мячей для гольфа из большой корзины с заменой. Размер и форма корзины зависят от предположения о распределении входных данных (например, нормальное распределение со средним значением 100 и стандартным отклонением 10, по сравнению с равномерным распределением или треугольным распределением), когда некоторые корзины глубже или более симметричны, чем другие, что позволяет вытаскивать определенные шары чаще, чем другие.

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

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

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

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

Понимание распределений вероятностей для моделирования методом Монте-Карло

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

Во-первых, вы собираете необработанные данные — в данном случае заработную плату каждого сотрудника отдела, не освобожденного от уплаты налогов. Во-вторых, вы организуете данные в осмысленный формат и отображаете данные в виде частотного распределения на диаграмме. Чтобы создать частотное распределение, разделите заработную плату на групповые интервалы и перечислите эти интервалы на горизонтальной оси диаграммы. Затем вы перечисляете количество или частоту сотрудников в каждом интервале на вертикальной оси диаграммы. Теперь вы можете легко увидеть распределение заработной платы, не освобожденной от уплаты налогов, внутри отдела.

Взглянув на график ниже, можно увидеть, что большинство сотрудников (примерно 60 из 180) зарабатывают от $7,00 до $9,00 в час.

Эти данные можно построить в виде распределения вероятностей. Распределение вероятностей показывает количество сотрудников в каждом интервале в виде доли от общего числа сотрудников. Чтобы создать распределение вероятностей, разделите количество сотрудников в каждом интервале на общее количество сотрудников и перечислите результаты на вертикальной оси диаграммы.

На приведенной ниже диаграмме показано количество сотрудников в каждой группе заработной платы в процентах от всех сотрудников; Можно оценить вероятность того, что сотрудник, выбранный случайным образом из всей группы, получит заработную плату в течение заданного интервала. Например, предполагая, что те же самые условия существовали на момент взятия выборки, вероятность того, что сотрудник, выбранный случайным образом из всей группы, составляет 0,33 (один шанс из трех), что он зарабатывает от 8,00 до 8,50 долларов в час.

Распределения вероятностей могут быть дискретными или непрерывными. Дискретные распределения вероятностей описывают различные значения, обычно целые числа, без промежуточных значений и отображаются в виде последовательности вертикальных столбцов. Дискретное распределение, например, может описывать количество орлов в четырех подбрасываниях монеты как 0, 1, 2, 3 или 4.

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

Выбор правильного распределения вероятностей

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

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

Равномерное распределение

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

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

Условия

Математические построения для равномерного распределения следующие:

  • Минимальное количество элементов фиксировано.
  • Максимальное количество элементов фиксировано.
  • Все значения между минимумом и максимумом встречаются с равной вероятностью.

Математические построения для равномерного распределения следующие:

Максимальное значение (Max) и минимальное значение (Min) являются параметрами распределения.

Требования к входным данным:

Min < Max и может принимать любое значение.

Предварительная обработка данных

Прежде всего, давайте импортируем все соответствующие библиотеки:import numpy as np
import pandas as pd
from scipy.stats import norm
from numpy.random import randn
from numpy import random as rn
import scipy.stats as si
from matplotlib import pyplot as plt
import seaborn as sns
from matplotlib import pyplot as plt
from IPython.display import Image
%matplotlib inline

Во-вторых, давайте создадим фрейм данных из 50 000 симуляций (т.е. испытаний) случайных чисел, которые равномерно распределены между 0 и 1.Simulations = 50000
f = np.arange(1,Simulations+1)
df = pd.DataFrame(data=f)
df[‘RAND’] = rn.rand(Simulations,1)
df.drop([0], axis=1, inplace=True)
df

И тогда для равномерного распределения есть формула

У нас есть общий диапазон и наша кумулятивная вероятность, которая на самом деле является нашими случайными числами (т.е. df["RAND"]). Итак, у нас есть генерация случайных чисел равномерного распределения.

Давайте просто отформатируем и это. Наш общий диапазон — это наш максимум минус наш минимум.TotalRange = Maximum — Minimum

Приведенная выше формула немного сложна для понимания, но давайте просто введем ее. Таким образом, нам нужно проверить, если кумулятивная вероятность (т.е. df[«df["RAND"][i]«][i]) очень близка к нулю (скажем, 0.001), то наша цена продажи (т.е. df["SalesPrice"][i] ) очень близка к нашему минимуму, если кумулятивная вероятность очень близка к единице (скажем, 0.999), то наша цена продажи очень близка к нашему максимуму.

И самый простой способ проверить, правильно ли мы написали нашу формулу, заключается в том, что если наша кумулятивная вероятность равна 0, нам нужно получить 8,00, что является минимумом, а если она равна 1, нам нужно получить 8,50, что является максимумом. Итак, наша формула верна.df[‘SalesPrice’] = 0.00
for i in range(0,len(df[‘SalesPrice’])):
df[‘SalesPrice’][i] = Minimum + df[‘RAND’][i]*(TotalRange)
df

Обозначим наш целевой столбец буквой V.V = df[‘SalesPrice’]
V

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

Входными данными функции является вероятность, представляющая область под стандартной нормальной кумулятивной функцией распределения (cdf) в диапазоне от минус бесконечности до плюс бесконечности. Выходными данными функции является число стандартных отклонений, соответствующих данной конкретной области при стандартной нормальной кумулятивной функции распределения.def NORMSINV(x):
x = si.norm.ppf(x)
return (x)

Например, обратное стандартному нормальному кумулятивному распределению, с вероятностью 0,9000 равно 1,2815515655446004NORMSINV(0.90)

1.2815515655446004

Анализ данных

Проанализируем информацию моделированияG = 857550884
np.random.seed(G)
from datetime import datetime
start_time = datetime.now()
from datetime import datetime
# datetime object containing current date and time
now = datetime.now()
dt_string = now.strftime(«%d/%m/%Y %H:%M:%S»)

time_elapsed1 = datetime.now() — start_time
from scipy import stats
Workbook_Name = «Uniform Distribution with MC Simulations.ipynb»
Number_of_Simulations = «{:,.0f}».format(Simulations)
Number_of_Iterations = «{:,.0f}».format(1)
Number_of_Inputs = «{:,.0f}».format(2)
Number_of_Outputs = 38
Sampling_Type = «Latin Hypercube»
Simulation_Start_Time = dt_string
Simulation_Duration = «{}».format(time_elapsed1)
Random_N_Generator = «Mersenne Twister»
Random_Seed = G
e = [«Workbook Name»,»Number of Simulations»,»Number of Iterations»,»Number of Inputs»,»Number of Outputs»,»Sampling Type»,\
«Simulation Start Time»,»Simulation Duration»,»Random # Generator»,»Random Seed»]
f = [Workbook_Name, Number_of_Simulations, Number_of_Iterations, Number_of_Inputs, Number_of_Outputs, Sampling_Type,\
Simulation_Start_Time, Simulation_Duration, Random_N_Generator, Random_Seed]
Per5 = «{:,.2f}».format(np.percentile(V, 5))
P5 = «{:.0%}».format(0.05)
Per10 = «{:,.2f}».format(np.percentile(V, 10))
P10 = «{:.0%}».format(0.10)
Per15 = «{:,.2f}».format(np.percentile(V, 15))
P15 = «{:.0%}».format(0.15)
Per20 = «{:,.2f}».format(np.percentile(V, 20))
P20 = «{:.0%}».format(0.20)
Per25 = «{:,.2f}».format(np.percentile(V, 25))
P25 = «{:.0%}».format(0.25)
Per30 = «{:,.2f}».format(np.percentile(V, 30))
P30 = «{:.0%}».format(0.30)
Per35 = «{:,.2f}».format(np.percentile(V, 35))
P35 = «{:.0%}».format(0.35)
Per40 = «{:,.2f}».format(np.percentile(V, 40))
P40 = «{:.0%}».format(0.40)
Per45 = «{:,.2f}».format(np.percentile(V, 45))
P45 = «{:.0%}».format(0.45)
Per50 = «{:,.2f}».format(np.percentile(V, 50))
P50 = «{:.0%}».format(0.50)
Per55 = «{:,.2f}».format(np.percentile(V, 55))
P55 = «{:.0%}».format(0.55)
Per60 = «{:,.2f}».format(np.percentile(V, 60))
P60 = «{:.0%}».format(0.60)
Per65 = «{:,.2f}».format(np.percentile(V, 65))
P65 = «{:.0%}».format(0.65)
Per70 = «{:,.2f}».format(np.percentile(V, 70))
P70 = «{:.0%}».format(0.70)
Per75 = «{:,.2f}».format(np.percentile(V, 75))
P75 = «{:.0%}».format(0.75)
Per80 = «{:,.2f}».format(np.percentile(V, 80))
P80 = «{:.0%}».format(0.80)
Per85 = «{:,.2f}».format(np.percentile(V, 85))
P85 = «{:.0%}».format(0.85)
Per90 = «{:,.2f}».format(np.percentile(V, 90))
P90 = «{:.0%}».format(0.90)
Per95 = «{:,.2f}».format(np.percentile(V, 95))
P95 = «{:.0%}».format(0.95)
Minimum = «{:,.2f}».format(np.min(V))
Maximum = «{:,.2f}».format(np.max(V))
Mean = «{:,.2f}».format(np.mean(V))
Std_Dev = «{:,.2f}».format(np.std(V))
Variance = np.var(V)
Std_Error = «{:,.4f}».format(np.std(V)/np.sqrt(Simulations))
Skewness = round(stats.skew(V),9)
Kurtosis = round((stats.kurtosis(V)+3),9)
Median = «{:,.2f}».format(np.median(V))
Mode = «{:,.2f}».format(stats.mode(V)[0][0])
Left_X = Per5
Left_P = P5
Right_X = Per95
Right_P = P95
Diff_X = «{:,.2f}».format((np.percentile(V, 95) — np.percentile(V, 5)))
Diff_P = «{:.0%}».format(0.90)
Confidence_Interval = P95
Lower_Bound = «{:,.2f}».format((np.mean(V) — (np.std(V)/np.sqrt(Simulations))*NORMSINV(0.975)))
Upper_Bound = «{:,.2f}».format((np.mean(V) + (np.std(V)/np.sqrt(Simulations))*NORMSINV(0.975)))
g = {«Information»: e, «Result»: f}
st = pd.DataFrame(data=g)
a = [«Minimum»,»Maximum»,»Mean»,»Std Dev»,»Variance»,»Std Error», «Skewness»,»Kurtosis»,»Median»,»Mode»,\
«Left X»,»Left P»,»Right X»,»Right P»,»Diff X»,»Diff P»,»Confidence Interval»,»Lower Bound»,»Upper Bound»]
b = [Minimum, Maximum, Mean, Std_Dev, Variance, Std_Error, Skewness, Kurtosis, Median, Mode, Left_X, Left_P,\
Right_X, Right_P, Diff_X, Diff_P, Confidence_Interval, Lower_Bound, Upper_Bound]
c = [P5,P10,P15,P20,P25,P30,P35,P40,P45,P50,P55,P60,P65,P70,P75,P80,P85,P90,P95]
d = [Per5, Per10, Per15, Per20, Per25, Per30, Per35, Per40, Per45, Per50, Per55, Per60, Per65,\
Per70, Per75, Per80, Per85, Per90, Per95]
d = {«Statistics»: a, «Statistics Result»: b, «Percentile»: c, «Percentile Result»: d}
st1 = pd.DataFrame(data=d)
from datetime import date
today = date.today()
now = datetime.now()
import calendar
curr_date = date.today()
print(«\033[1m Simulation Summary Information»)
print(«\033[0m ================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))
st

Проанализируем статистику по прогнозируемой цене продажиprint(«\033[1m Summary Statistics for Sales Price»)
print(«\033[0m ======================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))
st1

Визуализация данных

Визуализируем функцию плотности вероятности (pdf) для прогнозируемой цены продажиplt.figure(figsize = (4,4))
sns.set(font_scale = 1.2)
ax = sns.histplot(data=V,bins=10,color=»red»)
ax.set_xlabel(«Sales Price»,fontsize=14)
ax.set_ylabel(«Frequency»,fontsize=14)
ax.set_xlim([7.8,8.6])
print(«\033[1m Probability Density Function for Sales Price (Sim#1)»)
print(«\033[0m ======================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))

Визуализируем функцию кумулятивного распределения (cdf) для прогнозируемой цены продажиplt.figure(figsize = (4,4))
kwargs = {«cumulative»: True}
sns.axes_style(«whitegrid»)
sns.set(font_scale = 1.2)
ax = sns.distplot(V, hist_kws=None, kde_kws=kwargs,color=»red»)
ax.set_xlabel(«Sales Price»,fontsize=14)
ax.set_ylabel(«Frequency»,fontsize=14)
ax.set_xlim([7.8,8.6])
print(«\033[1m Cumulative Distribution Function for Sales Price (Sim#1)»)
print(«\033[0m ===========================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))

Заключение

Для того, чтобы предсказать цену продажи, я использовал моделирование по методу Монте-Карло и равномерное распределение. Моделирование методом Монте-Карло широко используется в ситуациях, когда отсутствует аналитическое решение и другие подходы (например, биномиальное дерево).

Моделирование по методу Монте-Карло использует моделирование факторов риска в соответствии с предполагаемым распределением вероятностей, а затем вычисляет цену продажи для каждого моделирования отдельно. Среднее значение цены продажи называется прогнозируемым значением цены продажи.

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

В прогнозах мы можем определить вероятность наступления, называемую доверительным интервалом. То есть, при двух значениях, каковы шансы, что результат окажется между этими двумя значениями?

В этих расчетах я использовал более 50 000 случайных симуляций. Существует 90% вероятность того, что окончательный результат (в данном случае цена продажи на 2023 год) будет между $8,03 и $8,47. С другой стороны, существует 5% вероятность того, что цена продажи будет ниже 8,03 доллара США, и еще 5% вероятность того, что цена продажи будет выше 8,47 доллара США.

То есть двусторонний доверительный интервал — это симметричный интервал, центрированный на среднем значении или медианном значении (т. е. значении 50-го процентиля). Таким образом, оба хвоста будут иметь одинаковую вероятность.

Среднее значение прогнозируемой цены продажи на 2023 год составляет 8,25 доллара США, что в точности является средним значением между минимальным значением прошлых данных о цене продажи (т. е. 8,00 долларов США) и максимальным значением прошлых данных о цене продажи (т. е. 8,25 доллара США).

Метод Монте-Карло: ретроспективное моделирование

Фото: Внутренняя стоимость

1. Мотивация

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

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

Упрощенно можно представить себе подход моделирования по методу Монте-Карло как многократное извлечение мячей для гольфа из большой корзины с заменой. Размер и форма корзины зависят от предположения о распределении входных данных (например, нормальное распределение со средним значением 100 и стандартным отклонением 10, по сравнению с равномерным распределением или треугольным распределением), когда некоторые корзины глубже или более симметричны, чем другие, что позволяет вытаскивать определенные шары чаще, чем другие.

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

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

В данной статье показано, как можно использовать ретроспективное моделирование Монте-Карло для прогнозирования неизвестного прошлого значения эквити с использованием логнормального распределения.

2. Марковское свойство

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

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

Статистические свойства истории цен на акции IBM могут быть полезны для определения характеристик стохастического процесса, за которым следует цена акции (например, ее волатильность). Здесь подчеркивается, что конкретный путь, по которому шли акции в прошлом, не имеет значения.

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

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

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

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

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

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

3. Стохастический процесс непрерывного времени

Рассмотрим переменную, которая следует марковскому стохастическому процессу. Предположим, что его текущее значение равно 10 и что изменение его значения в течение одного года равно φ(0, 1), где φ(μ, σ) обозначает распределение вероятностей, которое обычно распределяется со средним μ и стандартным отклонением σ.

Каково распределение вероятностей изменения значения переменной в течение двух лет? Изменение за два года является суммой двух нормальных распределений, каждое из которых имеет среднее значение, равное нулю, и стандартное отклонение, равное 1,0.

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

Дисперсия распределения вероятностей равна квадрату его стандартного отклонения. Таким образом, дисперсия изменения значения рассматриваемой нами переменной за один год равна 1,0

Таким образом, среднее значение изменения рассматриваемой нами переменной в течение двух лет равно нулю, а дисперсия этого изменения равна 2,0. Таким образом, изменение переменной за два года составляет φ(0, √2).

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

Мы предполагаем, что это одно и то же. Из этого следует, что дисперсия изменения в течение шестимесячного периода должна составлять 0,5. Эквивалентно, стандартное отклонение изменения равно √0,5, так что распределение вероятностей изменения значения переменной в течение шести месяцев равно φ(0, √0,5)

Аналогичный аргумент показывает, что изменение значения переменной в течение трех месяцев равно φ(0, √0,25). В более общем случае, изменение в течение любого периода времени длины T составляет φ(0, √T).

В частности, изменение длины δt за очень короткий период времени составляет φ(0, √δt). Знаки квадратного корня в этих результатах могут показаться странными.

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

Дисперсия изменения переменной в нашем примере равна 1,0 в год, так что дисперсия изменения за два года равна 2,0, а дисперсия изменения за три года равна 3,0. Стандартное отклонение изменения через два и три года составляет √2 и √3 соответственно.

Строго говоря, мы не должны говорить о стандартном отклонении переменной как о 1,0 в год. Это должно быть «1,0 на квадратный корень из лет». Полученные результаты объясняют, почему неопределенность часто называют пропорциональной квадратному корню из времени.

2.1. Винские процессы

Процесс, за которым следует рассматриваемая нами переменная, известен как винеровский процесс. Это особый тип марковского стохастического процесса со средним изменением, равным нулю, и коэффициентом дисперсии 1,0 в год.

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

Выражаясь формально, переменная z следует винеровскому процессу, если она обладает следующими двумя свойствами:

Свойство 1. Изменение δz за небольшой промежуток времени δt равно

где ε — случайный рисунок из стандартизированного нормального распределения, φ(0, 1).

Свойство 2. Значения δz для любых двух различных коротких интервалов времени δt независимы.

Из первого свойства следует, что δz само по себе имеет нормальное распределение с

Второе свойство подразумевает, что z следует марковскому процессу. Рассмотрим увеличение величины z в течение относительно длительного периода времени, T.

Это можно обозначить как z(T) — z(0). Его можно рассматривать как сумму приращений z за N малых интервалов времени длиныδ t, где

Таким образом

где ε(i = 1,2, . . . , N) — случайные рисунки из φ(0, 1). Из второго свойства винеровских процессов εi независимы друг от друга. Из приведенного выше уравнения следует, что z(T) — z(0) нормально распределяется с

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

По истечении пяти лет он обычно распределяется со средним значением 25 и стандартным отклонением √5, или 2,236. Наша неопределенность относительно значения переменной в определенный момент времени в будущем, измеряемая ее стандартным отклонением, возрастает как квадратный корень из того, как далеко мы смотрим вперед.

В обычном исчислении обычно исходят от небольших изменений к пределу по мере того, как малые изменения приближаются к нулю. Таким образом, δy/δx становится dy/dx в пределе, и так далее.

Аналогично можно поступить и при работе со стохастическими процессами. Винеровский процесс является пределом в виде δt → 0 описанного выше процесса для z.

На рисунке ниже показано, что происходит с траекторией, по которой следует z, при приближении к пределу δt → 0. Обратите внимание, что путь довольно «неровный».

Как получается винеровский процесс при δt → 0

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

  1. Ожидаемая длина пути, за которым следует z в любом интервале времени, бесконечна.
  2. Ожидаемое количество раз, когда z равно какому-либо конкретному значению в любом интервале времени, бесконечно.

2.2. Обобщенный винеровский процесс

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

Коэффициент дисперсии, равный 1,0, означает, что дисперсия изменения z во временном интервале длины T равна T. Обобщенный винеровский процесс для переменной x может быть определен в терминах dz следующим образом:

где a и b — константы.

Чтобы понять приведенное выше уравнение, полезно рассмотреть два компонента в правой части по отдельности. Член a dt подразумевает, что x имеет ожидаемую скорость дрейфа a в единицу времени.

Без члена b dz уравнение имеет вид

из чего следует, что

Интегрируя по времени, получаем

где x0 — значение x в нулевой момент времени. За период времени длины T величина x возрастает на величину aT.

Член b dz в правой части уравнения можно рассматривать как добавление шума или изменчивости к пути, по которому следует x. Величина этого шума или изменчивости равна b, умноженной на винеровский процесс.

Винский процесс имеет стандартное отклонение, равное 1,0. Из этого следует, что b, умноженное на винеровский процесс, имеет стандартное отклонение b.

На малом интервале времени δt изменение значения δx задается приведенными выше уравнениями как

где, как и раньше, ε — случайный отбор из стандартизированного нормального распределения. Таким образом, δx имеет нормальное распределение с

Аргументы, аналогичные аргументам, приведенным для винеровского процесса, показывают, что изменение значения x в любом интервале времени T нормально распределяется с

Таким образом, обобщенный винеровский процесс, приведенный в приведенном выше уравнении, имеет ожидаемый коэффициент дрейфа (т.е. средний дрейф в единицу времени) a и коэффициент дисперсии (т.е. дисперсию в единицу времени) . Это проиллюстрировано на рисунке ниже

Обобщенный винеровский процесс: a = 0,3, b = 1,5

Рассмотрим ситуацию, когда денежная позиция компании, измеряемая в тысячах долларов, подчиняется обобщенному процессу Винера с дрейфом 20 в год и коэффициентом отклонения 900 в год.

Изначально денежная позиция равна 50. В конце одного года денежная позиция будет иметь нормальное распределение со средним значением 70 и стандартным отклонением √900 или 30.

По истечении шести месяцев он будет иметь нормальное распределение со средним значением 60 и стандартным отклонением 30√0,5 = 21,21. Наша неопределенность в отношении денежной позиции в какой-то момент в будущем, измеряемая ее стандартным отклонением, возрастает по мере того, как квадратный корень из того, насколько далеко вперед мы смотрим.

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

2.3. Процесс Ито

Можно определить еще один тип стохастического процесса. Это называется процессом Ито. Это обобщенный винеровский процесс, в котором параметры a и b являются функциями величины базовой переменной x и времени t.

Алгебраически процесс Itô может быть записан

Как ожидаемая скорость дрейфа, так и степень дисперсии процесса Itô могут изменяться с течением времени. В малом промежутке времени между t и t + δt переменная изменяется от x до x + δx, где

Эта зависимость включает в себя небольшое приближение. Предполагается, что дрейф и дисперсия x остаются постоянными, равными a(x, t) и b(x, t)² соответственно, в течение интервала времени между t и t + δt.

3. Модель поведения цены акций

Модель поведения цены акций известна как геометрическое броуновское движение. Дискретно-временная версия модели

Переменная δS — это изменение цены акции S за небольшой интервал времени δt, а ε — случайная выборка из стандартизированного нормального распределения (т.е. нормального распределения со средним значением нуля и стандартным отклонением 1,0).

Параметр μ — это ожидаемая норма доходности в единицу времени от акции, а параметр σ — волатильность цены акции. Оба эти параметра считаются постоянными. Левая часть следующего уравнения

– доходность, обеспечиваемая запасом за короткий промежуток времени δt. Член μδt — это математическое ожидание этой доходности, а член σεδt — стохастическая составляющая доходности.

Дисперсия стохастической составляющей (и, следовательно, всей доходности) равна σ²δt. Это согласуется с определением волатильности как данности; То есть, σ такова, что σ√δt является стандартным отклонением доходности за короткий промежуток времени δt.

Приведенные выше уравнения показывают, как dS/S нормально распределяется со средним значением μ δt и стандартным отклонением σδt. Другими словами,

4. Свойство Lognormal

Теперь мы используем лемму Ито для получения процесса, за которым следует ln S, когда S следует следующему процессу

Определять

При постоянном μ и σ является разумной моделью движения цен на акции. Из леммы Ито следует, что процесс, за которым следует функция G из S и t, равен

Заметим, что и S, и G подвержены влиянию одного и того же основного источника неопределенности, dz. Это оказывается очень важным при выводе результатов Блэка-Шоулза.

Потому что

Из приведенного выше уравнения следует, что процесс, которому следует G, равен

Поскольку μ и σ постоянны, это уравнение показывает, что G = ln S подчиняется обобщенному винеровскому процессу. Имеет постоянную скорость дрейфа μ — σ²/2 и постоянную дисперсию σ².

Таким образом, изменение ln S между нулевым временем и некоторым будущим временем, T, нормально распределяется со средним значением

и дисперсия

Это означает, что

где ST — цена акции в будущий момент времени TS0 — цена акции в нулевой момент времени, а φ(m, s) обозначает нормальное распределение со средним m и стандартным отклонением s .

Следующее уравнение

показывает, что ln ST нормально распределена. Переменная имеет логнормальное распределение, если натуральный логарифм переменной нормально распределен.

Таким образом, модель поведения цены акций, которую мы разработали в этой статье, подразумевает, что цена акции в момент времени T, учитывая ее сегодняшнюю цену, логарифмически распределена. Стандартное отклонение логарифма цены акции составляет σt. Он пропорционален квадратному корню из того, как далеко вперед мы смотрим.

5. Реальный мир и нейтральный к риску мир

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

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

К сожалению, нелегко узнать правильную ставку дисконтирования, чтобы применить ее к ожидаемому выигрышу в реальном мире. Позиция по колл-опциону более рискованна, чем позиция по акциям, и предположим, что в реальном мире ожидаемая доходность по акциям составляет 16%.

В результате, ставка дисконтирования, применяемая к выплате по колл-опциону, превышает 16%. Не зная стоимости опциона, мы не знаем, насколько она должна быть больше 16%.

Риск-нейтральная оценка решает эту проблему. Мы знаем, что в мире с нейтральным риском ожидаемая доходность всех активов (и, следовательно, ставка дисконтирования, используемая для всех ожидаемых выплат) является безрисковой ставкой.

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

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

Для некоторого значения рыночной цены риска мы получаем «реальный мир» и наблюдаемые на практике темпы роста цен ценных бумаг.

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

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

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

Из этого следует, что актуарный подход верен только тогда, когда математическое ожидание от дериватива одинаково как в реальном, так и в нейтральном к риску мире.

Мы упоминали, что когда мы переходим из реального мира в мир нейтральный к риску, волатильность переменных остается прежней, но ожидаемые темпы их роста могут измениться. Например, ожидаемые темпы роста индекса фондового рынка уменьшаются, возможно, на 4% или 5%, когда мы переходим из реального мира в мир нейтральный к риску.

Кроме того, можно обоснованно предположить, что ожидаемый темп роста переменной одинаков как в реальном, так и в нейтральном к риску мире, если переменная имеет нулевой систематический риск (β=0), так что процентные изменения переменной имеют нулевую корреляцию с доходностью фондового рынка. Из этого можно сделать вывод, что актуарный подход к оценке производного инструмента дает правильный ответ, если все базовые переменные имеют нулевой систематический риск.

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

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

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

6. Перспективное моделирование по методу Монте-Карло

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

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

  1. Оцените ожидаемую доходность капитала в реальном мире для горизонта прогнозирования.
  2. Выборка случайного пути для S в реальном мире.
  3. Повторите шаг 2, чтобы получить множество примеров будущей стоимости акций в реальном мире.
  4. Вычислите среднее значение выборки будущей стоимости собственного капитала, чтобы получить оценку ожидаемой будущей стоимости собственного капитала.

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

где dz — это процесс Винера, μ — ожидаемая доходность в реальном мире, а σ — волатильность. Чтобы смоделировать путь, по которому следует S, мы можем разделить горизонт прогноза на N коротких интервалов длины Δt и получим

где S(t) обозначает значение S в момент времени tε — случайная выборка из нормального распределения со средним нулем и стандартным отклонением 1,0. Это позволяет вычислить значение S в момент времени Δ t из начального значения S, значение в момент времени 2Δ t вычислить из значения в момент времени Δ t и т. д.

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

Из леммы Ито процесс, за которым следует ln S, имеет вид

Таким образом, чтобы

или-эквивалентно

Это уравнение используется для построения пути для S.

Работа с ln S, а не с S дает большую точность. Кроме того, если μ и σ постоянны, то

справедливо для всех Т. Из этого следует, что

Это уравнение может быть использовано для прогнозирования будущей стоимости собственного капитала в момент времени T.

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

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

или

Поэтому я решил построить корректное уравнение для ретроспективного моделирования методом Монте-Карло с использованием логнормального распределения.

Ниже приведено уравнение Поланицера для прогнозирования прошлой стоимости капитала с использованием ретроспективного моделирования по методу Монте-Карло и логнормального распределения:

Где k решается итеративно так, что минимизирует разницу между оценкой прошлой стоимости справедливости в детерминированном мире (т.е. мире абсолютной достоверности) и оценкой прошлого значения справедливости в стохастическом мире (это среднее значение в момент времени 0 всех путей).

7. Пример ретроспективного моделирования методом Монте-Карло с использованием логнормального распределения

Предположим, что ожидаемая доходность собственного капитала (рассчитанная на основе CAPM) составляет 15% годовых (в реальном мире), а стандартное отклонение доходности (т.е. волатильность) составляет 50% годовых. Приведенная стоимость капитала в момент 0 составляет $211.700. Спрогнозируйте прошлую стоимость акций 5-летней давности и 95%-ный доверительный интервал для прошлой стоимости акций 5-летней давности.

Это означает, что ST = 211,700, μ = 0,15, σ = 0,50 и T = 5,00.

Во-первых, давайте загрузим соответствующие библиотеки pythonimport numpy as np
import pandas as pd
from scipy.stats import norm
from numpy.random import randn
from numpy import random as rn
from matplotlib import pyplot as plt
import seaborn as sns
%matplotlib inline

Во-вторых, давайте напишем наши входные переменныеST = 211.7000016612675
μ = 0.15
σ = 0.50
T = 5

В-третьих, давайте найдем прошлую стоимость акций 5-летней давности в условиях детерминированного мира (т.е. мира абсолютной определенности)

S0 = ST*np.exp(-μ*T)
S0

100.00000000000001

В-четвертых, решим использовать 50 000 путей и независимое от пути моделирование (т.е. N = 1, то есть разделим горизонт прогноза на N коротких интервалов длины Δt)Paths = 50000
N = 1
dt = T/N

В-пятых, давайте найдем прошлую стоимость акций 5-летней давности в стохастическом мире (т.е. мире неопределенности). Итак, мы начинаем со столбца ST, то есть с приведенной стоимости собственного капитала в момент времени 0data = ST*np.ones(Paths)
# initialize list of lists
data1 = {‘ST’:data}
# Create the pandas DataFrame
df = pd.DataFrame(data1)
# print dataframe.
df

Затем мы создадим столбец с именем RAND(), который генерирует случайные числа от 0 до 1.df[‘RAND()’] = rn.rand(Paths,1)
df

После этого мы создадим столбец с именем NORMSINV(), который дает обратную кумулятивную функцию для стандартного нормального распределения. Из этого следует, что NORMSINV(RAND()) дает случайную выборку из стандартного нормального распределения.df[‘NORMSINV()’] = norm.ppf(df[‘RAND()’])
df

Наконец, мы создадим столбец под названием S0то есть прошлую стоимость эквити в момент времени -T (5 лет назадЭто соответствует следующему уравнению и является случайной выборкой из множества для всех значений эквити:

Я итеративно искал и нашел, что k должно быть равно 1.5

df[‘S0’] = df[‘ST’]/np.exp(1.5-((μ-0.5*σ**2)*T+σ*df[‘NORMSINV()’]*np.sqrt(T)))
df

В-шестых, давайте проанализируем информацию моделированияV = df[‘S0’]

G = 1629562571
np.random.seed(G)
from datetime import datetime
start_time = datetime.now()
from datetime import datetime
# datetime object containing current date and time
now = datetime.now()
dt_string = now.strftime(«%d/%m/%Y %H:%M:%S»)
def NORMSINV(x):
x = si.norm.ppf(x)
return (x)
time_elapsed1 = datetime.now() — start_time
from scipy import stats
Workbook_Name = «Backward-Looking MC Simulation.ipynb»
Number_of_Simulations = «{:,.0f}».format(Paths)
Number_of_Iterations = «{:,.0f}».format(N)
Number_of_Inputs = «{:,.0f}».format(5)
Number_of_Outputs = 38
Sampling_Type = «Latin Hypercube»
Simulation_Start_Time = dt_string
Simulation_Duration = «{}».format(time_elapsed1)
Random_N_Generator = «Mersenne Twister»
Random_Seed = G
e = [«Workbook Name»,»Number of Simulations»,»Number of Iterations»,»Number of Inputs»,»Number of Outputs»,»Sampling Type»,\
«Simulation Start Time»,»Simulation Duration»,»Random # Generator»,»Random Seed»]
f = [Workbook_Name, Number_of_Simulations, Number_of_Iterations, Number_of_Inputs, Number_of_Outputs, Sampling_Type,\
Simulation_Start_Time, Simulation_Duration, Random_N_Generator, Random_Seed]
Per5 = «{:,.4f}».format(np.percentile(V, 5))
P5 = «{:.0%}».format(0.05)
Per10 = «{:,.4f}».format(np.percentile(V, 10))
P10 = «{:.0%}».format(0.10)
Per15 = «{:,.4f}».format(np.percentile(V, 15))
P15 = «{:.0%}».format(0.15)
Per20 = «{:,.4f}».format(np.percentile(V, 20))
P20 = «{:.0%}».format(0.20)
Per25 = «{:,.4f}».format(np.percentile(V, 25))
P25 = «{:.0%}».format(0.25)
Per30 = «{:,.4f}».format(np.percentile(V, 30))
P30 = «{:.0%}».format(0.30)
Per35 = «{:,.4f}».format(np.percentile(V, 35))
P35 = «{:.0%}».format(0.35)
Per40 = «{:,.4f}».format(np.percentile(V, 40))
P40 = «{:.0%}».format(0.40)
Per45 = «{:,.4f}».format(np.percentile(V, 45))
P45 = «{:.0%}».format(0.45)
Per50 = «{:,.4f}».format(np.percentile(V, 50))
P50 = «{:.0%}».format(0.50)
Per55 = «{:,.4f}».format(np.percentile(V, 55))
P55 = «{:.0%}».format(0.55)
Per60 = «{:,.4f}».format(np.percentile(V, 60))
P60 = «{:.0%}».format(0.60)
Per65 = «{:,.4f}».format(np.percentile(V, 65))
P65 = «{:.0%}».format(0.65)
Per70 = «{:,.4f}».format(np.percentile(V, 70))
P70 = «{:.0%}».format(0.70)
Per75 = «{:,.4f}».format(np.percentile(V, 75))
P75 = «{:.0%}».format(0.75)
Per80 = «{:,.4f}».format(np.percentile(V, 80))
P80 = «{:.0%}».format(0.80)
Per85 = «{:,.4f}».format(np.percentile(V, 85))
P85 = «{:.0%}».format(0.85)
Per90 = «{:,.4f}».format(np.percentile(V, 90))
P90 = «{:.0%}».format(0.90)
Per95 = «{:,.4f}».format(np.percentile(V, 95))
P95 = «{:.0%}».format(0.95)
Minimum = «{:,.4f}».format(np.min(V))
Maximum = «{:,.4f}».format(np.max(V))
Range = «{:,.4f}».format((np.max(V))-(np.min(V)))
Count = len(V)
Mean = «{:,.4f}».format(np.mean(V))
Median = «{:,.4f}».format(np.median(V))
Mode = «{:,.4f}».format(stats.mode(V)[0][0])
Variance = np.var(V)
Std_Dev = «{:,.4f}».format(np.std(V))
Std_Err = «{:,.4f}».format(np.std(V)/np.sqrt(Count))
Skewness = round(stats.skew(V),9)
Kurtosis = round((stats.kurtosis(V)+3),9)
Exc_Kurt = round((stats.kurtosis(V)),9)
Left_X = Per5
Left_P = P5
Right_X = Per95
Right_P = P95
Diff_X = «{:,.4f}».format((np.percentile(V, 95) — np.percentile(V, 5)))
Diff_P = «{:.0%}».format(0.90)
g = {«Information»: e, «Result»: f}
st = pd.DataFrame(data=g)
a = [«Minimum»,»Maximum»,»Range»,»Count»,»Mean»,»Median»,»Mode»,»Variance»,»Std Dev»,»Std Err»,»Skewness»,»Kurtosis»\
,»Exc_Kurt»,»Left X»,»Left P»,»Right X»,»Right P»,»Diff X»,»Diff P»]
b = [Minimum, Maximum, Range, Count, Mean, Median, Mode, Variance, Std_Dev, Std_Err, Skewness, Kurtosis\
,Exc_Kurt, Left_X, Left_P, Right_X ,Right_P, Diff_X, Diff_P]
c = [P5,P10,P15,P20,P25,P30,P35,P40,P45,P50,P55,P60,P65,P70,P75,P80,P85,P90,P95]
d = [Per5, Per10, Per15, Per20, Per25, Per30, Per35, Per40, Per45, Per50, Per55, Per60, Per65,\
Per70, Per75, Per80, Per85, Per90, Per95]
d = {«Statistics»: a, «Statistics Result»: b, «Percentile»: c, «Percentile Result»: d}
st1 = pd.DataFrame(data=d)
from datetime import date
today = date.today()
now = datetime.now()
import calendar
curr_date = date.today()
print(«\033[1m Simulation Summary Information»)
print(«\033[0m ================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))
st

В-седьмых, проанализируем статистику по прошлому значению эквитиprint(«\033[1m Summary Statistics for the Past Value of Equity»)
print(«\033[0m ================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))
st1

В-восьмых, давайте визуализируем функцию плотности вероятности (pdf) для прошлого значения эквитиplt.figure(figsize = (4,4))
sns.set(font_scale = 1.2)
ax = sns.histplot(data=V,bins=500,color=»red»)
ax.set_xlabel(«Values»,fontsize=14)
ax.set_ylabel(«Frequency»,fontsize=14)
ax.set_xlim(np.percentile(V, 2.5),np.percentile(V, 97.5))
print(«\033[1m Probability Density Function for the Past Value of Equity (Sim#1)»)
print(«\033[0m =================================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))

Наконец, давайте визуализируем функцию кумулятивного распределения (cdf) для прошлой стоимости собственного капиталаplt.figure(figsize = (4,4))
kwargs = {«cumulative»: True}
sns.axes_style(«whitegrid»)
sns.set(font_scale = 1.2)
ax = sns.distplot(V, hist_kws=None, kde_kws=kwargs,color=»red»)
ax.set_xlabel(«Values»,fontsize=14)
ax.set_ylabel(«Frequency»,fontsize=14)
ax.set_xlim(np.percentile(V, 2.5),np.percentile(V, 97.5))
print(«\033[1m Cumulative Distribution Function for the Past Value of Equity (Sim#1)»)
print(«\033[0m =====================================================================»)
print(«\033[1m Performed By:»,»\033[0mRoi Polanitzer»)
print(«\033[1m Date:»,»\033[0m»,calendar.day_name[curr_date.weekday()],»,»,today.strftime(«%B %d, %Y»),»,»,now.strftime(«%H:%M:%S AM»))

Теперь найдем оценку будущей стоимости собственного капиталаV = df[‘ST’]
print(«The past value of equity is: {:.2f}».format(np.mean(V)))

Прошлая стоимость собственного капитала: 99.96

Давайте оценим точность нашей оценкиprint(«The stochastic error of the past value of equity is: {:.2f}».format(np.std(V)/np.sqrt(Paths)))

Стохастическая ошибка прошлого значения эквити составляет: 0,71

Стандартные отклонения нормы, соответствующие 95% доверительному интервалу, составляют αMIN = −1,96 и αMAX = 1,96. Другими словами, у нас есть 2,5% в каждом хвосте.print(«The lower bound for the past value of equity with 95% confidence level is {:.4}».format(np.mean(V) — 1.96*np.std(V)/np.sqrt(Paths)))

Нижняя граница для прошлой стоимости акций с доверительной вероятностью 95% составляет 98,5664print(«The upper bound for the past value of equity with 95% confidence level is {:.4}».format(np.mean(V) + 1.96*np.std(V)/np.sqrt(Paths)))

Верхняя граница для прошлой стоимости акций с доверительной вероятностью 95% составляет 101,3594

8 Заключение

Мы видим, что наша стохастическая мировая оценка прошлой стоимости собственного капитала ($99,9629) очень близка к нашей детерминированной мировой оценке прошлой стоимости капитала ($100,00). Это означает, что мы строим хорошую модель для прогнозирования прошлых значений капитала, используя ретроспективное моделирование по методу Монте-Карло.

Ценообразование опционов по методу Монте-Карло

Найти теоретическую цену опционов с помощью моделирования по методу Монте-Карло

  • Сложность:★★★★☆
  • Варианты цен с помощью моделирования по методу Монте-Карло
  • Внедрить методику снижения отклонений для улучшения результата прайнинга

Предисловие

Моделирование методом Монте-Карло получило широкое применение в области финансовых исследований. В разделе 【Quant(19)】Prediction of Portfolio Performance мы рассказали, как использовать моделирование по методу Монте-Карло для прогнозирования цены акций. В сегодняшней статье мы расширим применение на более сложное ценообразование опционов. В 【Quant】CRR Model и 【Quant】Модель Блэка-Шоулза и греки мы использовали соответственно принципы биномиальных деревьев и формулу Блэка-Шоулза для расчета теоретических цен опционов. В этих статьях также разъяснялись многие фундаментальные концепции, связанные с опционами. Для тех, кто имеет ограниченное представление о вариантах, рекомендуется сначала прочитать эти две статьи, прежде чем переходить к текущей. В последующих разделах мы продемонстрируем, как использовать моделирование по методу Монте-Карло для прогнозирования цен на акции. Затем мы расширим дискуссию до ценообразования европейских опционов и, наконец, представим некоторые методы снижения дисперсии, которые помогут нам использовать моделирование по методу Монте-Карло.

Среда программирования и требуемый модуль

В качестве редактора используется Windows 11 и Jupyter Notebook# Import required packages
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm
import time
import tejapi
plt.style.use(‘bmh’)

# Log in TEJ API
api_key = ‘YOUR_KEY’
# tejapi.ApiConfig.api_base=»http://10.10.10.66″
tejapi.ApiConfig.api_key = api_key
tejapi.ApiConfig.ignoretz = True

База данных

База данных по торговле акциями: Нескорректированная ежедневная цена акций, код базы данных (TWN/APRCD).
База данных деривативов: ежедневная информация о сделках опционов, код базы данных (TWN/AOPTION).

Импорт данных

С использованием нескорректированных цен закрытия Тайваньского взвешенного фондового индекса (Y9999) за период с 31 января 2021 года по 19 апреля 2023 года. Мы также загрузим опционы колл и пут взвешенного индекса Тайваня (TXO202304C15500, TXO202304P15500). Это опционы европейского типа, с датой начала торгов 31 января и датой экспирации 19 апреля. Цена исполнения установлена на уровне 15 500. Кроме того, мы установим столбец «mdate» (дата) в качестве индекса.# Import required data
gte, lte = ‘2021-03-16’, ‘2023-04-20’
stocks = tejapi.get(‘TWN/APRCD’, # stock price
paginate = True,
coid = ‘Y9999’,
mdate = {‘gte’:gte, ‘lte’:lte},
opts = {
‘columns’:[ ‘mdate’,’close_d’]
}
)
# Get options price
puts = tejapi.get( # puts price
‘TWN/AOPTION’,
paginate = True,
coid = ‘TXO202304P15500’,
mdate = {‘gte’:gte, ‘lte’:lte},
opts = {
‘columns’:[‘mdate’, ‘coid’,’settle’, ‘kk’, ‘theoremp’, ‘acls’, ‘ex_price’, ‘td1y’, ‘avolt’, ‘rtime’]
}
)
calls = tejapi.get( # calls price
‘TWN/AOPTION’,
paginate = True,
coid = ‘TXO202304C15500’,
mdate = {‘gte’:gte, ‘lte’:lte},
opts = {
‘columns’:[‘mdate’, ‘coid’,’settle’, ‘kk’, ‘theoremp’, ‘acls’, ‘ex_price’, ‘td1y’, ‘avolt’, ‘rtime’]
}
)

# set index
stocks = stocks.set_index(‘mdate’)
puts = puts.set_index(‘mdate’)
calls = calls.set_index(‘mdate’)

Обработка данных

Рассчитывая дневную доходность и перемещая волатильность, установите в качестве окна 252 дня.# Calculate daily return
stocks[‘daily return’] = np.log(stocks[‘close_d’]) — np.log(stocks[‘close_d’].shift(1))
stocks[‘moving volatility’] = stocks[‘daily return’].rolling(252).std()

Прогноз цены акций

Мы используем моделирование по методу Монте-Карло для моделирования траекторий цен на акции. Концепция моделирования по методу Монте-Карло довольно проста. Он включает в себя получение процесса возврата актива и его дискретизацию, а затем использование небольших временных интервалов для расчета изменений цен на активы. Например, рассматривая цены на акции, их доходность следует геометрическому броуновскому движению. Таким образом, мы можем получить дискретизированное стохастическое дифференциальное уравнение (уравнение 1), где Wt представляет собой винеровский процесс. Применив формулу Ито, мы получаем уравнение 2 в качестве основного уравнения для моделирования по методу Монте-Карло для прогнозирования цен на акции, где Zt следует стандартному нормальному распределению.

Далее мы можем использовать Python для программирования приведенного выше уравнения. Суть метода Монте-Карло заключается в одновременной оценке нескольких траекторий цены акций с помощью уравнения 2. Наконец, суммируя и усредняя последнюю цену акции по каждому пути, мы получаем прогнозируемую цену акции. Здесь мы определяем следующие переменные:

  • S0: Цена актива на начальную дату
  • r: Историческая доходность актива
  • sigma: стандартное отклонение исторической доходности актива
  • T: Время созревания
  • Nsteps: Номера временных интервалов
  • Nrep: Номера стоковых путей

Кодифицируем приведенное выше уравнение в функцию «mc_asset» для выполнения моделирования по методу Монте-Карло.def mc_asset(S0, r, sigma, T, Nsteps, Nrep):
SPATH = np.zeros((Nrep, 1 + Nsteps))
SPATH[:, 0] = S0
dt = T / Nsteps
nudt = (r — 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)

for i in range(0,Nrep):
for j in range(0,Nsteps):
SPATH[i,j+1] = SPATH[i,j] * np.exp(nudt + sidt * np.random.normal())
return SPATH

После настройки функции мы можем установить аргументы, свидетельствующие о том, может ли функция работать или нет. Результат можно представить на рисунке 1. Каждая строка на рисунке 1 представляет собой траекторию моделирования запаса.S0 = 100
K = 110
CallOrPut = ‘call’
r = 0.03
sigma = 0.25
T = 0.5
Nsteps = 10000
Nrep = 1000
SPATH = mc_asset(S0, r, sigma, T, Nsteps, Nrep)

plt.figure(figsize = (10,8))
for i in range(len(SPATH)):
plt.plot(SPATH[i])
plt.xlabel(‘Numbers of steps’)
plt.ylabel(‘Stock price’)
plt.title(‘Monte Carlo Simulation for Stock Price’)
plt.show()

Рисунок 1: Моделирование по методу Монте-Карло на цене акций

Ценообразование опционов

Мы можем использовать вышеупомянутый метод для прогнозирования цены акций на момент экспирации опционов. Затем мы вычисляем внутреннюю стоимость опциона по истечении срока действия для каждого пути. Наконец, дисконтировав внутреннюю стоимость обратно к текущей стоимости и взяв среднее значение, мы можем получить теоретическую цену опциона. Пожалуйста, обратитесь к конкретному коду ниже:def mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep):
SPATH = mc_asset(S0, r, sigma, T, Nsteps, Nrep)
if CallOrPut == ‘call’:
payoffs = np.maximum(SPATH[:,-1] — K, 0)
return np.mean(payoffs)*np.exp(-r*T)
else:
payoffs = np.maximum(K — SPATH[:,-1], 0)
return np.mean(payoffs)np.exp(-rT)

Мы можем рассчитать цену колла или пута, просто изменив аргумент «CallOrPut». Далее мы предоставляем всем аргументам соответствующее значение для проверки. Более того, мы сравниваем результат с теоретической ценой, рассчитанной на основе 【Quant】 модели Блэка Шоулза и греков.S0 = 100
K = 110
CallOrPut = ‘put’
r = 0.03
sigma = 0.25
T = 0.5
Nsteps = 10000
Nrep = 1000
p_ = mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep)
mybs = BS_formula(S0, K, r, sigma, T)
c, p = mybs.BS_price()

print(f’Monte Carlo price: {c_}’)
print(f’Black Scholes price: {p}’)

Результаты показаны ниже на рисунке 2. Как видите, разница между ценой моделирования по методу Монте-Карло и формулой Блэка-Шоулза невелика.

Рисунок 2: Сравнение цены МС и БС

Антитезный вариативный метод

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

Концепция, лежащая в основе этого метода, заключается в создании траектории цены для базового актива (уравнение 3) и одновременной генерации траектории с противоположной доходностью (уравнение 4). В этом случае корреляция между двумя путями равна -1, что приводит к минимальной ковариации при объединении двух путей. Это, в свою очередь, снижает волатильность при оценке цены опциона.

Конкретная реализация кода выглядит следующим образом. Мы генерируем две матрицы для вычислений: SPATH1 представляет прямой путь, а SPATH2 представляет путь в обратном направлении. Можно заметить, что в каждой итерации при генерации случайных чисел (эпсилон) обе матрицы имеют одинаковые случайные числа. В результате уменьшаются вычислительные усилия, а также снижается волатильность при прогнозировании цен опционов.def mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep):

SPATH1 = np.zeros((int(Nrep/2), 1 + Nsteps))
SPATH2 = np.zeros((int(Nrep/2), 1 + Nsteps))
SPATH1[:, 0], SPATH2[:, 0] = S0, S0
dt = T / Nsteps
nudt = (r — 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)

for i in range(0,int(Nrep/2)):
for j in range(0,Nsteps):
epsilon = np.random.normal()
SPATH1[i,j+1] = SPATH1[i,j] * np.exp(nudt + sidt * epsilon)
SPATH2[i,j+1] = SPATH2[i,j] * np.exp(nudt — sidt * epsilon)

if CallOrPut == ‘call’:
C1 = np.maximum(SPATH1[:, -1] — K, 0)
C2 = np.maximum(SPATH2[:, -1] — K, 0)
C = np.mean(0.5 * (C1 + C2))
C0 = np.exp(-r*T) * C
return C0
else:
P1 = np.maximum(K — SPATH1[:, -1], 0)
P2 = np.maximum(K — SPATH2[:, -1], 0)
P = np.mean(0.5 * (P1 + P2))
P0 = np.exp(-r*T) * P
return P0

Затем мы вводим значения и проверяем, согласуются ли цены опционов, полученные с помощью метода антитетической переменной, с результатами обычного моделирования по методу Монте-Карло. Результаты можно увидеть на рисунке 3, и можно заметить, что они действительно очень близки друг к другу.CallOrPut = ‘put’
K = 110
S0 = 100
r = 0.03
sigma = 0.25
T = 0.5
Nrep = 10000
Nsteps = 1000
print(‘Price under AV: ‘, mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print(‘Price under MC: ‘, mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))

Рисунок 3: Сравнение AV и MC

Метод управления переменной

В дополнение к методу Antithetic Variate, мы также можем уменьшить волатильность теоретических цен опционов с помощью метода Control Variate. Предположим, что у нас есть две случайные величины, X и Y, где вычисление среднего значения и дисперсии переменной Y является простым. Предположим, что эти две переменные можно объединить, чтобы сформировать новую переменную Z (уравнение 5). В этом случае математическое ожидание Z совпадает с математическим ожиданием X (как показано в уравнении 6), а дисперсия определяется параметром c. Таким образом, мы можем найти оптимальное значение для c*, которое минимизирует дисперсию Z, как показано в уравнении 7. Мы рассматриваем X и Y как опционные и базовые цены акций для каждого пути. Используя ковариацию между историческими ценами опционов и акций и дисперсию цен акций, мы можем вычислить оптимальное c*, а затем вычислить теоретическую цену опциона (E(Z)), используя уравнение 5. Это позволяет достичь цели снижения волатильности в симуляциях по методу Монте-Карло.

Результат программирования отображается следующим образом:

  • CallOrPut: выбор между коллом и путом
  • K: Цена исполнения
  • S0: Цена актива в начальный момент времени
  • r: Историческая доходность актива
  • sigma: стандартное отклонение исторической доходности актива
  • T: Время созревания
  • Nsteps: Номера временных интервалов
  • Nrep: Номера стоковых путей
  • Npilot: Период времени для получения оптимального c*

def mc_options_CV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep, Npilot):

# Calculate covariance between stock and options price
SPATH = np.zeros((Npilot, 1 + Nsteps))
SPATH[:, 0] = S0
dt = T / Nsteps
nudt = (r — 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)

for i in range(0,Npilot):
for j in range(0,Nsteps):
SPATH[i,j+1] = SPATH[i,j] * np.exp(nudt + sidt * np.random.normal())
Sn = SPATH[:, -1]
if CallOrPut == ‘call’:
Cn = np.maximum(SPATH[:,-1] — K, 0) * np.exp(-r*T)
MatCov = np.cov(Sn, Cn)[0,1]
VarY = S0 ** 2 * np.exp(2 * r * T) * (np.exp(T * sigma ** 2) — 1)
c = -MatCov / VarY
ExpY = S0 * np.exp(r*T)
else:
Pn = np.maximum(K — SPATH[:,-1], 0) * np.exp(-r*T)
MatCov = np.cov(Sn, Pn)[0,1]
VarY = S0 ** 2 * np.exp(2 * r * T) * (np.exp(T * sigma ** 2) — 1)
c = -MatCov / VarY
ExpY = S0 * np.exp(r*T)

# Applying control variate function with optimal c*
SPATH2 = np.zeros((Nrep, 1 + Nsteps))
SPATH2[:, 0] =S0
dt = T / Nsteps
nudt = (r — 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)

for i in range(0,Nrep):
for j in range(0,Nsteps):
SPATH2[i,j+1] = SPATH2[i,j] * np.exp(nudt + sidt * np.random.normal())
S = SPATH2[:, -1]
if CallOrPut == ‘call’:
C = np.maximum(SPATH2[:,-1] — K, 0) * np.exp(-r*T)
CVC = np.mean(C + c * (S — ExpY))
return CVC
else:
P = np.maximum(K — SPATH2[:,-1], 0) * np.exp(-r*T)
CVP = np.mean(P + c * (S — ExpY))
return CVP

Вводя числа в аргументы, мы свидетельствуем о разнице между ценой опционов из метода Antithetic Variate и Control Variate. Результат показан на рисунке 4. Можно сказать, что они идентичны.

Рисунок 4: Сравнение AV и CV

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

Таким образом, мы изучили три метода расчета теоретических цен опционов с помощью моделирования по методу Монте-Карло: традиционный метод, метод антитетической переменной и метод контрольной вариации. Далее мы представим реальные примеры для сравнения теоретических коэффициентов, полученных с помощью каждого метода, с коэффициентами Блэка-Шоулза, рассчитанными TEJ, чтобы увидеть, есть ли между ними существенные различия.S0 = stocks.loc[‘2023-01-31’][‘close_d’]
K = 15500
r = stocks[‘daily return’].rolling(252).mean().loc[‘2023-01-31’] # average return of stock
T = 51 / 252
sigma = stocks.loc[‘2023-01-31’][‘moving volatility’] * np.sqrt(252)
Nstep = 50000
Nrep = 50000
Npilot = 5000
CallOrPut = ‘put’

print(‘Monte Carlo theoretical price: ‘, mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print(‘Monte Carlo with AV theoretical price: ‘, mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print(‘Monte Carlo with CV theoretical price: ‘, mc_options_CV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep, Npilot))
print(‘TEJ Black Scholes price: ‘, puts.loc[‘2023-01-31’][‘theoremp’])
print(‘Real price: ‘, puts.loc[‘2023-01-31’][‘settle’])

Мы используем пут-опцион Тайваньской фондовой биржи с ценой исполнения 15 500 и периодом времени с 1 января 2023 года по 19 апреля 2023 года. Сигма рассчитывается как стандартное отклонение доходности Тайваньской фондовой биржи с окном в 252 дня, а r — средняя доходность за последние 252 дня. Мы предполагаем, что сегодня 31 января 2023 года. Результаты показаны на рисунке 5, где мы видим, что цены, полученные с помощью трех методов, ближе к фактическим ценам по сравнению с рассчитанными ценами TEJ.

Рисунок 5: Сравнение для всех методов и цены TEJ BS

Заключение

Метод ценообразования по методу Монте-Карло в большей степени опирается на закон больших чисел по сравнению с моделью CRR и моделью Блэка-Шоулза. Он постепенно аппроксимирует разумные теоретические цены, моделируя большое количество траекторий цен на акции. Ожидается, что в сегодняшнюю эпоху значительного прогресса в производительности компьютеров алгоритмы и модели ценообразования, основанные на данных, будут становиться все более распространенными. Занимаясь торговлей опционами, инвесторы также могут рассмотреть возможность включения метода Монте-Карло в свои соображения. В будущем мы предоставим больше знаний, связанных с опционами и производными продуктами. Пожалуйста, продолжайте следить за нами. Кроме того, читатели и инвесторы могут приобрести решения от TEJ E Shop для создания собственных программ ценообразования опционов. Обратите внимание, что упомянутые выше формулы ценообразования и опционные продукты предназначены только для демонстрационных целей и не являются нашими инвестициями или рекомендациями.

Исходный код

https://medium.com/media/37469de7c5960cd0ee396a28f2f5f3b8

Расширенное чтение

Ссылки по теме

Игровая площадка квантов

Представьте, что вы пытаетесь предсказать, сколько раз монета упадет орлом после того, как подбросит ее 10 раз. Вы можете попробовать сложную математику, но более простой способ — подбросить монету 10 раз и посмотреть, что произойдет!

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

В этой статье мы раскроем секреты этой увлекательной техники Монте-Карло, изучим ее работу, применение в квантовых финансах и то, какую пользу она может принести вам как кванту.

Понимание математики, стоящей за волшебством:

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

  • Нормальное распределение: Для переменных с симметричным колоколообразным распределением, таких как ежедневные изменения цен.
  • Логнормальное распределение: Для переменных с положительным перекосом, таких как цены на активы, на более длительных горизонтах.
  • Биномиальное распределение: Для событий с двумя возможными исходами, таких как успех/неудача сделки.

Используя эти распределения, мы создаем «пути», представляющие различные потенциальные эволюции финансовых переменных. Каждый путь действует как виртуальный сценарий, разворачивающийся в течение дискретных временных шагов. Ключевые уравнения включают:

  • Динамика цен: S(t+1) = S(t) * exp((r — d²/2) * dt + σ * ε * sqrt(dt)), где S — цена, r — безрисковая ставка, d — дивидендная доходность, σ — волатильность, ε — случайное число из выбранного распределения, dt — шаг по времени.
  • Ценообразование опционов: Формула Блэка-Шоулза или другие модели ценообразования опционов используются для расчета цен опционов на основе смоделированных цен на активы.

Зачем использовать моделирование по методу Монте-Карло?

Есть несколько веских причин, по которым моделирование по методу Монте-Карло является основным элементом инструментария Quant:

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

Применение в Quant Finance:

Универсальность моделирования по методу Монте-Карло обеспечивает широкий спектр применений в квантовых финансах:

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

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

Представьте себе, что менеджер хедж-фонда оценивает сложную опционную стратегию, включающую несколько базовых активов. Используя моделирование по методу Монте-Карло, они могут:

  1. Определите модель с помощью уравнений, представляющих динамику ценкорреляцию между активами и доходность опционов.
  2. Присваивайте распределения вероятностей соответствующим переменным, таким как волатильность и процентные ставки.
  3. Запускайте тысячи симуляций, генерируя различные рыночные сценарии на основе случайных чисел.
  4. Анализируйте результаты, рассчитывая ожидаемую доходность, потенциальные убытки и статистику рисков, такую как Value at Risk (VaR).

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

Реализация стандартного варианта использования на Python:

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

Метод:

  1. Моделирование траекторий движения акций по ценам:
  • Разделяйте время на маленькие шаги.
  • Генерируйте случайные числа, представляющие ежедневные изменения цен, используя геометрическую формулу броуновского движения, которая включает в себя волатильность и случайность.
  • Моделируйте тысячи различных ценовых траекторий для акции в течение срока действия опциона.

2. Рассчитать выплату опциона:

  • Для каждого смоделированного пути рассчитайте выплату опциона при наступлении срока погашения на основе окончательной цены акции и цены исполнения. Это просто максимум разницы между ценой акции и ценой исполнения (ноль, если он отрицательный).

3. Ориентировочная цена опциона:

  • Усредняйте выигрыши опционов по всем смоделированным траекториям.
  • Дисконтируйте эту среднюю выплату обратно к настоящему, используя безрисковую ставку и время до погашения. Это дает ориентировочную цену колл-опциона.
import numpy as np
import matplotlib.pyplot as plt

# Define option parameters
S0 = 100 # Initial stock price
K = 105 # Strike price
T = 1 # Time to maturity
sigma = 0.2 # Volatility
r = 0.05 # Risk-free rate

# Define number of simulations
n_sims = 10000

# Generate random price paths using geometric Brownian motion
dt = T / 252 # Discretize time into daily steps
dW = np.random.normal(loc=0, scale=np.sqrt(dt), size=(n_sims, 252)) # Generate random increments
S = S0 * np.exp((r - 0.5 * sigma**2) * dt + sigma * dW).cumsum(axis=1)

# Calculate option payoffs at maturity
call_payoff = np.maximum(S[:, -1] - K, 0)

# Calculate option price using Monte Carlo average
call_price = np.exp(-r * T) * np.mean(call_payoff)

# Print option price
print("Call option price:", call_price)

# Plot distribution of simulated call payoffs
plt.hist(call_payoff, bins=50, density=True, label="Call Payoff Distribution")
plt.xlabel("Call Payoff")
plt.ylabel("Density")
plt.title("Distribution of Simulated Call Payoffs")
plt.grid(True)
plt.show()

Этот код python создает следующий граф:

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

Имитационное моделирование методом Монте-Карло

Знакомство

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

Шаг 1: Настройка среды

Для начала нам нужно установить необходимые библиотеки Python. К таким pandasnumpyhttpxbacktestingpandas_tamatplotlibscipyrich и другие. Вот как их установить:

import pandas as pd
import numpy as np
from datetime import datetime
import concurrent.futures
import warnings
from rich.progress import track
from backtesting import Backtest, Strategy
import pandas_ta as ta
import matplotlib.pyplot as plt
from scipy.stats import norm
import httpx

warnings.filterwarnings("ignore")

Шаг 3: Определение служебных функций

Нам нужны функции для получения исторических цен на акции и криптовалюты из API:

def make_api_request(api_endpoint, params):
with httpx.Client() as client:
# Make the GET request to the API
response = client.get(api_endpoint, params=params)
if response.status_code == 200:
return response.json()
print("Error: Failed to retrieve data from API")
return None

def get_historical_price_full_crypto(symbol):
api_endpoint = f"{BASE_URL_FMP}/historical-price-full/crypto/{symbol}"
params = {"apikey": FMP_API_KEY}
return make_api_request(api_endpoint, params)


def get_historical_price_full_stock(symbol):
api_endpoint = f"{BASE_URL_FMP}/historical-price-full/{symbol}"
params = {"apikey": FMP_API_KEY}

return make_api_request(api_endpoint, params)

def get_SP500():
api_endpoint = "https://en.wikipedia.org/wiki/List_of_S%26P_500_companies"
data = pd.read_html(api_endpoint)
return list(data[0]['Symbol'])

def get_all_crypto():
"""
All possible crypto symbols
"""
return [
"BTCUSD",
"ETHUSD",
"LTCUSD",
"BCHUSD",
"XRPUSD",
"EOSUSD",
"XLMUSD",
"TRXUSD",
"ETCUSD",
"DASHUSD",
"ZECUSD",
"XTZUSD",
"XMRUSD",
"ADAUSD",
"NEOUSD",
"XEMUSD",
"VETUSD",
"DOGEUSD",
"OMGUSD",
"ZRXUSD",
"BATUSD"
"USDTUSD",
"LINKUSD",
"BTTUSD",
"BNBUSD",
"ONTUSD",
"QTUMUSD",
"ALGOUSD",
"ZILUSD",
"ICXUSD",
"KNCUSD",
"ZENUSD",
"THETAUSD",
"IOSTUSD",
"ATOMUSD",
"MKRUSD",
"COMPUSD",
"YFIUSD",
"SUSHIUSD",
"SNXUSD",
"UMAUSD",
"BALUSD",
"AAVEUSD",
"UNIUSD",
"RENBTCUSD",
"RENUSD",
"CRVUSD",
"SXPUSD",
"KSMUSD",
"OXTUSD",
"DGBUSD",
"LRCUSD",
"WAVESUSD",
"NMRUSD",
"STORJUSD",
"KAVAUSD",
"RLCUSD",
"BANDUSD",
"SCUSD",
"ENJUSD",
]

def get_financial_statements_lists():
api_endpoint = f"{BASE_URL_FMP}/financial-statement-symbol-lists"
params = {"apikey": FMP_API_KEY}
return make_api_request(api_endpoint, params)

Шаг 4: Разделение данных на обучающий и тестовый наборы

Далее мы получим исторические цены акций по данному символу и разобьем данные на обучающую и тестовую выборки:

stock_symbol = "AAPL"
stock_prices = get_historical_price_full_stock(stock_symbol)
data = pd.DataFrame(stock_prices['historical'])

# Splitting the data
january_2023_index = data[(data['date'] < '2023-01-01')].index
prices_after_january_2023 = data.drop(january_2023_index)

# Assuming prices_after_january_2023 is your DataFrame
prices_after_january_2023 = prices_after_january_2023.rename(columns={
'open': 'Open',
'high': 'High',
'low': 'Low',
'close': 'Close',
'volume': 'Volume' # Only if you have a volume column
})

# Drop any additional columns that are not required
required_columns = ['date', 'Open', 'High', 'Low', 'Close', 'Volume']
prices_after_january_2023 = prices_after_january_2023[required_columns]

# sort by date ascending
prices_after_january_2023 = prices_after_january_2023.sort_values(by=['date'], ascending=True)
plt.figure(figsize=(10, 6))
plt.title('Stock Prices')
plt.xlabel('Date')
plt.ylabel('Price')
plt.plot(prices_before_january_2023['date'], prices_before_january_2023['Close'])
plt.plot(prices_after_january_2023['date'], prices_after_january_2023['Close'])
plt.show()

Шаг 5: Моделирование по методу Монте-Карло

Мы выполним моделирование по методу Монте-Карло для прогнозирования будущих цен на акции:

def monte_carlo_simulation(data, days, iterations):
# Convert to numpy array if data is a pandas Series
if isinstance(data, pd.Series):
data = data.to_numpy()

# Ensure data is a numpy array
if not isinstance(data, np.ndarray):
raise TypeError("Data must be a numpy array or pandas Series")

# Calculate log returns
log_returns = np.log(data[1:] / data[:-1])
mean = np.mean(log_returns)
variance = np.var(log_returns)
drift = mean - (0.5 * variance)
daily_volatility = np.std(log_returns)

future_prices = np.zeros((days, iterations))
current_price = data[-1]
for t in range(days):
shocks = drift + daily_volatility * norm.ppf(np.random.rand(iterations))
future_prices[t] = current_price * np.exp(shocks)
current_price = future_prices[t]
return future_prices

Визуализация

simulation_days = 364
mc_iterations = 1000
start_state = strata.iloc[-1]
mc_prices = monte_carlo_simulation(prices_before_january_2023['Close'], simulation_days, mc_iterations)

# Last closing price repeated for each iteration
last_close_price = np.full((1, mc_iterations), prices_before_january_2023['Close'].iloc[-1])

# Combine the last closing price with the Monte Carlo simulation prices
mc_prices_combined = np.concatenate((last_close_price, mc_prices), axis=0)
markov_states = markov_chain_simulation(transition_matrix, start_state, simulation_days)

# Adjust the periods in the date range to match the number of rows in mc_prices_combined
last_date = prices_before_january_2023['date'].iloc[-1]
simulated_dates = pd.date_range(start=last_date, periods=simulation_days + 1)

# Visualizing the Monte Carlo simulation alongside historical data
plt.figure(figsize=(10, 6))

# Plot historical data
plt.plot(pd.to_datetime(prices_before_january_2023['date']), prices_before_january_2023['Close'], label='Before January 2023', linewidth=2)
plt.plot(pd.to_datetime(prices_after_january_2023['date']), prices_after_january_2023['Close'], label='After January 2023', linewidth=2)

# Taking average of all simulations on the 365th day
future_price_mcmc = np.mean(mc_prices_combined[364, :])
print(f"Average future price after 364 days: ${future_price_mcmc:.2f}")
print("Date of simulation:", simulated_dates[364].date())
simulated_date = simulated_dates[364].date()
real_price = prices_after_january_2023['Close'][18]
real_date = prices_after_january_2023['date'][18]
print(f"Real price on {real_date}: ${real_price:.2f}")
print(f"Real date: {real_date}")
print(f"Price difference: ${future_price_mcmc - real_price:.2f}")
print(f"Price difference percentage: {(future_price_mcmc - real_price) / real_price * 100:.2f}%")
print(f"Model accuracy: {100 - abs((future_price_mcmc - real_price) / real_price * 100):.2f}%")
# Plot Monte Carlo simulations mean
plt.plot(simulated_dates, mc_prices_combined.mean(axis=1), label='Monte Carlo Mean', linewidth=3)


# Plot each simulation path
for i in range(mc_iterations):
plt.plot(simulated_dates, mc_prices_combined[:, i], linewidth=0.5, color='gray', alpha=0.01)

plt.title('Monte Carlo Simulation of Stock Prices')
plt.xlabel('Date')
plt.ylabel('Price')
plt.legend()
plt.show()

Заключение

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

https://github.com/AlgoETS/MarkokChainMonteCarlo/tree/main

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник

Источник