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

"Начиная"

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

Могут ли продажи ванильного мороженого обогнать шоколад?

Содержание:

Вступление
Постановка задачи
Подготовка данных
Неправильный метод 1 - Независимое моделирование (параметрическое)
Неправильный метод 2 - Независимое моделирование (непараметрическое)
Метод 1 - многомерное распределение
Метод 2 - Копулы с маргинальными распределениями
Метод 3 - моделирование исторических комбинаций роста продаж
Метод 4— Декорреляция роста продаж магазина с помощью PCA

Вступление

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

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

В этой статье мы рассмотрим проблему коррелированных переменных в моделировании Монте-Карло. Мы рассмотрим 4 подходящих подхода для обработки корреляции.

Постановка задачи

В нашем выборочном наборе данных 20 лет назад годовые продажи шоколадного мороженого составляли 5 миллионов долларов, а продажи ванили - 4 миллиона долларов. Спрос на шоколадное мороженое был на 25% больше, чем на ванильное.

Спрос на шоколад и ваниль вырос за предыдущие 20 лет, но разрыв между продажами шоколада и ванили сократился. В прошлом году продажи шоколадного мороженого составили 12,6 миллиона долларов, ванили - 11,9 миллиона долларов. Спрос на шоколад сейчас только на 6% превышает спрос на ваниль.

ПРИМЕЧАНИЕ. Продажи шоколада и ванильного мороженого положительно коррелируют друг с другом, т. е. продажи шоколадного мороженого и ванильного мороженого имеют тенденцию увеличиваться вместе и уменьшаться вместе.

Вопрос. Может ли следующий год стать годом, когда продажи ванильного мороженого превысят продажи шоколадного мороженого?

Решение: мы можем запустить моделирование Монте-Карло, чтобы ответить на наш вопрос.

Было бы НЕПРАВИЛЬНО моделировать продажи шоколадного и ванильного мороженого независимо друг от друга.

Подготовка данных

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

START_CHOC = 5
START_VAN = 4

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

Мы предполагаем, что средняя скорость роста для обоих вкусов составляет 5% со стандартным отклонением 10%, а корреляция между скоростью роста обоих вкусов составляет 0,9.

# Config for multivariate normal distribution
CORR = 0.9
MU = 1.05
SIGMA = 0.1

Ниже приведен код для расчета производимых темпов роста и годовых продаж.

# Generate sales and growth rates for 20 periods
# using a multivariate normal distribution
cov_matrix = [[(SIGMA**2), CORR*(SIGMA**2)], 
              [CORR*(SIGMA**2), (SIGMA**2)]]
distribution = ss.multivariate_normal(
    mean=[MU, MU], cov=cov_matrix)
sample_data = distribution.rvs(20, random_state=5)
chocolate_growth = sample_data[:, 0]
vanilla_growth = sample_data[:, 1]
chocolate_sales = np.cumprod([START_CHOC] + list(chocolate_growth))
vanilla_sales = np.cumprod([START_VAN] + list(vanilla_growth))
# Prepare dataframes
growth_df = pd.DataFrame(data={
    "Growth Rate - Chocolate": chocolate_growth,
    "Growth Rate - Vanilla": vanilla_growth})
sales_df = pd.DataFrame(data={
    "Sales Chocolate (USD mln)": chocolate_sales,
    "Sales Vanilla (USD mln)": vanilla_sales})
df = pd.concat([growth_df, sales_df], axis=1)
df

Неправильный метод 1 - Независимое моделирование (параметрическое)

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

Мы знаем параметры распределения, поскольку мы создали набор данных, предполагая, что темпы роста продаж как шоколада, так и ванили следует нормальному распределению со средним значением 1,05 (рост в среднем 5%) и стандартным отклонением 0,1 (10%). См. Предыдущий раздел.

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

# Estimate means
mean_choc = np.mean(df["Growth Rate - Chocolate"])
print(f"Mean (Chocolate): {mean_choc}")
mean_van = np.mean(df["Growth Rate - Vanilla"])
print(f"Mean (Vanilla): {mean_van}")
# Estimate standard deviations
std_choc = np.std(df["Growth Rate - Chocolate"],
                 ddof=1)
print(f"StDev (Chocolate): {std_choc}")
std_van = np.std(df["Growth Rate - Vanilla"],
                ddof=1)
print(f"StDev (Vanilla): {std_van}")
# Define normal distributions
distribution_choc = ss.norm(mean_choc, std_choc)
distribution_van = ss.norm(mean_van, std_van)

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

growth_vanilla = distribution_choc.rvs(1000, 
    random_state=1)
growth_chocolate = distribution_van.rvs(1000, 
    random_state=2)

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

def exceed_check(growth_vanilla, growth_chocolate):
    '''This function takes sample growth rates for
    vanilla and chocolate ice cream sales and checks
    if vanilla ice cream sales would exceed chocolate
    given the combination of growth rates.
    
    Args:
    growth_vanilla (list): growth rates for vanilla sales
    growth_chocolate (list): growth rates for chocolate sales
    
    Returns:
    flags_list (list): A list of True/False flags indicating
        whether vanilla sales exceeds chocolate sales for the
        given combination of growth rates
    mean_pass_rate: Mean value of flags_list indicating the 
        probability vanilla sales exceeds chocolates given the 
        combination of growth rates
    '''
    # Last year's chocolate and vanilla ice cream sales
    # set as constants to improve performance
    FINAL_CHOCOLATE = 12.59
    FINAL_VANILLA = 11.94
    
    flags_list = [v*FINAL_VANILLA > c*FINAL_CHOCOLATE 
        for v,c in zip(growth_vanilla, growth_chocolate)]
    
    mean_pass_rate = np.mean(flags_list)
    
    return flags_list, mean_pass_rate

Наконец, давайте проверим вероятность того, что продажи ванили превысят шоколад…

exceed_check(growth_vanilla, growth_chocolate)[1]

Вероятность: 34,1%

Неправильный метод 2 - Независимое моделирование (непараметрическое)

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

growth_vanilla = df["Growth Rate - Vanilla"].dropna(
    ).sample(1000, replace=True, random_state=1)
growth_chocolate = df["Growth Rate - Chocolate"].dropna(
    ).sample(1000, replace=True, random_state=2)

Теперь мы проверяем вероятность превышения продаж ванили над шоколадом…

exceed_check(growth_vanilla, growth_chocolate)[1]

Вероятность: 35,9%

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

Метод 1 - многомерное распределение

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

Помните, что мы фактически сгенерировали данные из многомерного нормального распределения со средним значением 1,05 как для шоколада, так и для ванили, стандартными отклонениями 0,1 и корреляцией 0,9. Однако здесь мы предположим, что не знаем параметров, и будем использовать исторические данные для оценки параметров распределения, как мы это делали в независимом параметрическом моделировании (раздел «Неправильный метод 1»).

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

corr = df.corr().iloc[0, 1]
corr

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

# We need to recalculate the covariance matrix
# using the estimated paramaters
cov_matrix = [[std_choc**2, corr*std_choc*std_van], 
              [corr*std_choc*std_van, std_van**2]]
growth_rates = ss.multivariate_normal(
    mean=[mean_choc, mean_van],
    cov=cov_matrix).rvs(1000, random_state=1)
growth_chocolate = growth_rates[:, 0]
growth_vanilla = growth_rates[:, 1]

Наконец, мы вычисляем вероятность превышения продаж ванили над шоколадом…

exceed_check(growth_vanilla, growth_chocolate)[1]

Вероятность: 11,3%! Обратите внимание, как вероятность теперь намного ниже после того, как мы приняли во внимание корреляцию между продажами шоколада и ванильного мороженого.

Похоже, шоколад все-таки может остаться доминирующим ароматом ...

Метод 2 - Копулы с маргинальными распределениями

В предыдущем разделе мы рассмотрели многомерное нормальное распределение. Продажи как шоколада, так и ванильного мороженого следовали нормальному распределению, и они коррелированы. Это упростило представление их в многомерном нормальном распределении.

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

Подробное объяснение связок выходит за рамки, но проще говоря ...

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

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

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

corr = ss.kendalltau(
    df.dropna().iloc[:, 0], df.dropna().iloc[:, 1])[0]
corr

Теперь мы можем построить нашу связку ...

# We set variances to 1 because the covariance matrix we 
# are constructing will be used with a multivariate normal 
# distribution of means 0 and std 1 to derive a copula
cov_matrix = [[1, corr], 
              [corr, 1]]
# We will draw 1000 combinations from the distribution
random_vals = ss.multivariate_normal(cov=cov).rvs(
    1000, random_state=1)
# Finally a cumulative density function for a distribution
# ranges between 0 to 1 so it will be used to convert the
# drawn samples to the values we need for our copula
copula = ss.norm.cdf(random_vals)
print(copula.shape)
copula

Мы видим, что связка имеет 1000 коррелированных записей по 2 переменным, и они варьируются от 0 до 1.

Ниже приведен график разброса коррелированных значений копулы для двух наших переменных.

sns.scatterplot(x=copula[:, 0], y=copula[:, 1])

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

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

# distribution_choc and distribution_van area already
# calculated in previous sections
growth_chocolate = distribution_choc.ppf(copula[:, 0])
growth_vanilla = distribution_van.ppf(copula[:, 1])

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

exceed_check(growth_vanilla, growth_chocolate)[1]

Вероятность: 11,1%

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

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

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

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

Метод 3 - Моделирование исторических комбинаций роста продаж

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

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

Мы воспроизводим историю и берем темпы роста для обоих вкусов из 5-го года в первом раунде моделирования, например, затем за 11-й год во 2-м раунде и т. Д.

В независимом подходе в разделе «Неправильный метод 2» мы могли бы взять темп роста для ванили, например, с 5-го года, и скорость роста для шоколада с 11-го года в том же раунде моделирования. Мы переигрывали историю продаж шоколада и ванильного мороженого независимо друг от друга.

resampled_history = df[["Growth Rate - Chocolate", 
    "Growth Rate - Vanilla"]].dropna(
    ).sample(1000, replace=True, random_state=1)
growth_chocolate = resampled_history["Growth Rate - Chocolate"]
growth_vanilla = resampled_history["Growth Rate - Vanilla"]
exceed_check(growth_vanilla, growth_chocolate)[1]

Вероятность: 21%. Эта вероятность выше ожидаемых ~ 11%, полученных из двух последних подходящих методов, но ниже ~ 35%, предполагаемых неправильными подходами к моделированию.

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

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

Метод 4 - Декорреляция роста продаж магазина с помощью PCA

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

В независимом непараметрическом моделировании (раздел «Неправильный метод 2») подход оказался неверным, поскольку мы рассматривали темпы роста продаж ванили и шоколада как независимые. Это позволило нам получить 400 (20 x 20) точек исторических данных для повторной выборки, поскольку мы могли комбинировать любой из 20 исторических темпов роста продаж шоколада с любым из 20 исторических темпов роста ванили. Мы смогли совместить темпы роста шоколада в 1-м году с темпами роста ванили во 2-м году… в методе 3 мы не смогли.

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

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

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

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

pca = PCA(n_components=2)
# Normalize chocolate and vanilla growth rates to apply PCA
normalized_chocolate = (
    df["Growth Rate - Chocolate"].dropna() - mean_choc
) / std_choc
normalized_vanilla = (
    df["Growth Rate - Vanilla"].dropna() - mean_van
) / std_van
# Apply PCA transformation on normalized values
components = pca.fit_transform(
    np.array([normalized_chocolate, 
    normalized_vanilla]).T)
components

Затем мы запускаем 1000 симуляций, в которых отбираем образцы из массива компонентов.

np.random.seed(1)
sampled_components = [[x[0], x[1]] for x in zip(
    # Sampling 1000 entries from first PCA component
    np.random.choice(components[:, 0], 1000),
    # 2nd PCA component
    np.random.choice(components[:, 1], 1000))]

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

inverse_pca = pca.inverse_transform(sampled_components)
# Denormalizing
growth_chocolate = [(x * std_choc) + mean_choc 
    for x in inverse_pca[:, 0]]
growth_vanilla = [(x * std_van) + mean_van 
    for x in inverse_pca[:, 1]]

Наконец-то…

exceed_check(growth_vanilla, growth_chocolate)[1]

Вероятность: 17,1%, что ближе к 11%, которые мы ожидаем на основе методов 1 и 2.

Заключение

Ниже приведены вероятности, которые мы создали для того, чтобы продажи ванильного мороженого превысили объем продаж шоколада. Мы использовали 2 неправильных подхода и 4 подходящих.

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

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

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

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

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

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

"Начиная"