Распределение вероятностей количества дефолтов портфеля: итеративный подход

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

 В ней он, в частности, утверждает, что

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

Нужно понимать, что автор —  не специалист в области финансов. Но Владимир движется в правильном направлении. Он старается копать все глубже, а в результате может получить неожиданный и полезный результат. Его глаз не замылен академическими формулировками. У него нет цели рекламировать и монетизировать свой тг-канал, как это делают уже примелькавшиеся на смартлабе финансовые блогеры. Он хочет честно во всем разобраться с точки зрения обычного инвестора, не желающего терять лишние деньги на рынке ВДО. Он чем-то напоминает Алексея Галицкого, который, имея бэкграунд школьного учителя ОБЖ,  благодаря упорному труду создал свое «народное» рейтинговое агентство и стал заметной фигурой в области оценки кредитных рисков. Однако пока что к изысканиям Владимира стоит относиться с осторожностью: его анализ в текущем виде — это ещё не экспертное заключение. В общем, будем наблюдать. 

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

Допустим, у нас в портфеле имеется n облигаций и заданы вероятности дефолтов на определенном горизонте (cumulative default probabilities) 

⠀⠀⠀⠀⠀⠀⠀⠀q₁,  q₂ ,…, qₙ. 

Мы хотим найти P(L = k)   — вероятность того, что за время удержания портфеля произойдет ровно k дефолтов. 

Если предположить, что события дефолтов независимы, то случайная величина L (общее число дефолтов) подчиняется так называемому биномиальному распределению Пуассона (Poisson binomial distribution). Это распределение суммы независимых, но не одинаково распределённых бернуллиевских случайных величин. Мы можем построить его итеративно, добавляя по одной облигации, последовательно вычисляя распределение вероятностей числа дефолтов.

Обозначим:

⠀⠀⠀⠀⠀⠀pₖ[ j ] = P(произошло k дефолтов среди первых j облигаций) 

Вероятность иметь  k дефолтов после добавления j-й облигации равна сумме вероятностей следующих  событий

•  среди первых j − 1 облигаций было k дефолтов и в новой облигации не случился дефолт. 

• среди первых  j − 1 облигаций было k − 1 дефолтов и в новой  облигации случился дефолт:

⠀⠀⠀⠀⠀⠀⠀⠀pₖ[ j ]=pₖ[ j−1]⋅(1−qⱼ) + pₖ₋₁[ j−1]⋅qⱼ

Начальные условия (для первой облигации): 

                 p₀[1] = 1- q₁,       p₁[1] = q₁

После прохода по всем n облигациям получаем полный вектор распределения вероятностей pₖ,   k = 0, …, n.

Из полученного распределения легко найти необходимые метрики портфеля: вероятности того, что потери превысят заданный порог, квантили для расчета Value-at-Risk, условный VaR (Expected Shortfall) и т.д. Заметим, что нам не нужно знать вид распределения  P(L = k) для определения ожидаемых потерь портфеля в силу линейности математического ожидания. Для этого достаточно исходного набора вероятностей q₁, q₂,…,qₙ.

Ниже можно скачать файл Excel для расчета вероятностей.

Ссылка на файл  Excel на Яндекс.диск.

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

Мы можем учесть корреляцию дефолтов через однофакторную модель гауссовской копулы. В этой модели вводится общий фактор  —  нормально распределенная случайная величина  M ∼ N(0,1),  а условная вероятность дефолта j-й облигации при заданном значении M  выражается через ее безусловную вероятность qⱼ   и параметр корреляции ρ:

 Тогда общее распределение количества дефолтов получается усреднением условных распределений по фактору M.

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm


# --- INPUT ---
q = np.array([0.1,0.09,0.06,0.12,0.14,0.15,
              0.2,0.03,0.02,0.1,0.11,0.09,
              0.08,0.1,0.1,0.02,0.17,0.15,0.1,0.08])

N = len(q)
rho = 0.2   # корреляция
A = 1       # номинал
R = 0.4     # recovery


# --- Poisson binomial ---
def poisson_binomial(q):
    N = len(q)
    p = np.zeros(N+1)
    p[0] = 1 - q[0]
    p[1] = q[0]

    for i in range(1, N):
        new_p = np.zeros(N+1)
        for k in range(i+2):
            if k == 0:
                new_p[k] = p[k] * (1 - q[i])
            else:
                new_p[k] = p[k]*(1 - q[i]) + p[k-1]*q[i]
        p = new_p

    return p


# --- conditional default prob ---
def q_conditional(q, rho, M):
    return norm.cdf((norm.ppf(q) - np.sqrt(rho)*M) / np.sqrt(1 - rho))


# --- integrate over M ---
M_grid = np.linspace(-5, 5, 2000)
dM = M_grid[1] - M_grid[0]

p_total = np.zeros(N+1)

for M in M_grid:
    qM = q_conditional(q, rho, M)
    p_cond = poisson_binomial(qM)
    p_total += p_cond * norm.pdf(M) * dM


# --- normalize ---
p_total /= np.sum(p_total)


# --- metrics ---
k = np.arange(N+1)

expected_defaults = np.sum(k * p_total)
variance_defaults = np.sum((k - expected_defaults)**2 * p_total)

loss_per_default = A * (1 - R)
expected_loss = expected_defaults * loss_per_default
variance_loss = variance_defaults * loss_per_default**2


# --- VaR and CVaR ---
alpha = 0.95  # confidence level

# For number of defaults
var_k = np.quantile(k, alpha, method='linear')
cvar_k = np.mean(k[k >= var_k])

# For losses in monetary terms
losses = k * loss_per_default
var_loss = np.quantile(losses, alpha, method='linear')

# CVaR with probability weighting (more accurate)
tail_indices = losses >= var_loss
tail_probs = p_total[tail_indices]
if np.sum(tail_probs) > 0:
    cvar_loss = np.sum(losses[tail_indices] * tail_probs) / np.sum(tail_probs)
else:
    cvar_loss = var_loss


print("Expected defaults:", expected_defaults)
print("Variance defaults:", variance_defaults)
print("Expected loss:", expected_loss)
print("Variance loss:", variance_loss)

print(f"\n--- Risk metrics (α = {alpha:.0%}) ---")
print(f"VaR (defaults): {var_k:.2f} defaults")
print(f"CVaR (defaults): {cvar_k:.2f} defaults")
print(f"VaR (loss): {var_loss:.4f}")
print(f"CVaR (loss): {cvar_loss:.4f}")


# --- plots ---
plt.figure()
plt.bar(k, p_total)
plt.xlabel("Number of defaults")
plt.ylabel("Probability")
plt.title("Default distribution")
plt.show()


# cumulative loss distribution
cdf = np.cumsum(p_total)

plt.figure()
plt.plot(k, cdf)
plt.xlabel("Number of defaults")
plt.ylabel("CDF")
plt.title("Cumulative distribution")
plt.show()

Владимир в своей статье исследует насущную проблему  — как получить распределение доходности портфеля  ВДО-облигаций? Он делает достаточно сильные предположения:

1. Все облигации удерживаются до погашения.

2. Купоны реинвестируются в ту же самую облигацию, которая их выплатила.

3. Доходность к погашению каждой облигации постоянна. 

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

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

Поэтому подход Владимира применим именно к портфелю бескупонных облигаций. В этом случае:

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

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

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

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

import matplotlib.pyplot as plt
from scipy.stats import norm
from collections import defaultdict

# --- INPUT ---
y = np.array([0.2608,0.2295,0.2640,0.2926,0.2599,
              0.2839,0.2434,0.2608,0.3531,0.2799])

q = np.array([0.0855,0.0704,0.0855,0.0704,0.0704,
              0.1249,0.1036,0.0855,0.1249,0.1790])

w = np.ones(len(y)) / len(y)

T = 3
rho = 0.0  # можно менять



FV = (1 + y)**T

# точность дискретизации
STEP = 0.005


def round_key(x):
    return round(x / STEP) * STEP


# --- DP распределение при фиксированном q ---
def portfolio_distribution(q_vec):
    dist = {0.0: 1.0}

    for i in range(len(q_vec)):
        new_dist = defaultdict(float)
        payoff = w[i] * FV[i]

        for v, p in dist.items():
            # default
            new_dist[round_key(v)] += p * q_vec[i]

            # no default
            new_dist[round_key(v + payoff)] += p * (1 - q_vec[i])

        dist = dict(new_dist)

    return dist


# --- интегрирование по фактору ---
M_grid = np.linspace(-5, 5, 120)
dM = M_grid[1] - M_grid[0]

final_dist = defaultdict(float)

for M in M_grid:
    # условные вероятности дефолта
    qM = norm.cdf((norm.ppf(q) - np.sqrt(rho)*M) / np.sqrt(1 - rho))

    dist_M = portfolio_distribution(qM)

    weight = norm.pdf(M) * dM

    for v, p in dist_M.items():
        final_dist[v] += p * weight


# --- нормировка ---
total = sum(final_dist.values())
for k in final_dist:
    final_dist[k] /= total


# --- перевод в массивы ---
V_vals = np.array(list(final_dist.keys()))
P_vals = np.array(list(final_dist.values()))

# --- HPR ---
returns = V_vals**(1/T) - 1


# --- АГРЕГАЦИЯ ---
agg = defaultdict(float)

for r, p in zip(returns, P_vals):
    key = round(r, 4)
    agg[key] += p

returns_agg = np.array(list(agg.keys()))
P_agg = np.array(list(agg.values()))

# --- сортировка ---
idx = np.argsort(returns_agg)
returns_agg = returns_agg[idx]
P_agg = P_agg[idx]


# --- метрики ---
E_V = np.sum(V_vals * P_vals)
Var_V = np.sum((V_vals - E_V)**2 * P_vals)

E_R = np.sum(returns_agg * P_agg)
Var_R = np.sum((returns_agg - E_R)**2 * P_agg)
# --- CDF ---
cdf = np.cumsum(P_agg)


# --- VaR и Expected Shortfall ---
alpha = 0.95

# VaR
VaR_idx = np.searchsorted(cdf, 1 - alpha)
VaR = returns_agg[VaR_idx]

# Expected Shortfall
ES = np.sum(returns_agg[:VaR_idx] * P_agg[:VaR_idx]) / np.sum(P_agg[:VaR_idx])

print("Expected return:", E_R)
print("Variance return:", Var_R)
print(f"VaR (95%): {VaR}")
print(f"Expected Shortfall (95%): {ES}")

Expected return: 0.226337870388432
Variance return: 0.002048024847409808

VaR (95%): 0.1319
Expected Shortfall (95%): 0.10605073884031972

# Задаем границы диапазонов как переменные

LOWER_BOUND_1 = 0.14
UPPER_BOUND_1 = 0.18

LOWER_BOUND_2 = 0.20
UPPER_BOUND_2 = 0.27

TARGET_BOUND = 0.27

# 1. Средневзвешенная YTM (доходность без учета дефолтов)
expected_ytm = np.sum(y * w)

# 2. Вероятность потери капитала (доходность < 0)
neg_prob = np.sum(P_agg[returns_agg < 0])

# 3. Вероятность доходности < 14%
low_prob = np.sum(P_agg[returns_agg < LOWER_BOUND_1])

# 4. Вероятность доходности в диапазоне
mid_low_prob = np.sum(P_agg[(returns_agg >= LOWER_BOUND_1) & (returns_agg < UPPER_BOUND_1)])

# 5. Вероятность доходности в диапазоне
mid_prob = np.sum(P_agg[(returns_agg >= LOWER_BOUND_2) & (returns_agg < UPPER_BOUND_2)])

# 6. Вероятность доходности ≥ целевой
high_prob = np.sum(P_agg[returns_agg >= TARGET_BOUND])

# --- ВЫВОД С АВТОМАТИЧЕСКИМ ФОРМАТИРОВАНИЕМ ---
print("\n" + "="*60)
print("Анализ распределения доходности портфеля (точный расчёт)")
print("="*60)

print(f"Ожидаемая доходность без учёта риска (средневзвешенная YTM): {expected_ytm:.2%}")
print(f"Ожидаемая доходность с учётом дефолтов: {E_R:.2%}")


print(f"\n1. Потери капитала (доходность < 0): {neg_prob:.4%}")
print(f"   → {neg_prob:.2%} инвесторов теряют часть или все вложенные средства.")

print(f"\n2. Доходность < {LOWER_BOUND_1:.0%}: {low_prob:.4%}")
print(f"\n3. Доходность в диапазоне {LOWER_BOUND_1:.0%}–{UPPER_BOUND_1:.0%}: {mid_low_prob:.4%}")
print(f"\n4. Доходность в диапазоне {LOWER_BOUND_2:.0%}–{UPPER_BOUND_2:.0%}: {mid_prob:.4%}")
print(f"\n5. Целевая доходность (≥ {TARGET_BOUND:.0%}): {high_prob:.4%}")
print("\n" + "="*60)

============================================================
Анализ распределения доходности портфеля (точный расчёт)
============================================================
Ожидаемая доходность без учёта риска (средневзвешенная YTM): 27.28%
Ожидаемая доходность с учётом дефолтов: 22.63%

1. Потери капитала (доходность < 0): 0.0394%
   → 0.04% инвесторов теряют часть или все вложенные средства.

2. Доходность < 14%: 6.7848%

3. Доходность в диапазоне 14%–18%: 6.9944%

4. Доходность в диапазоне 20%–27%: 39.0139%

5. Целевая доходность (≥ 27%): 34.6258%

============================================================