Пенсионные рецепты на языке R

Прочитал тут на досуге книгу Moshe A. Milevsky «Retirement Income Recipes in R». В расчете познакомиться с языком R и понять как актуарная наука (страховщики) работает со случайностью времени жизни, на чем автор довольно подробно останавливался, а я раньше игнорировал. И как это можно включать в симуляции. В этом (весьма техническом) посте хочу поделиться некоторыми полезностями оттуда. Может, кому-то пригодятся.

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

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

\theta(t)=\Pr\bigl[T_x\ge t\bigr]=\exp\!\Big\{e^{\frac{x-m}{b}}\bigl(1-e^{\frac{t}{b}}\bigr)\Big\},

где $T_x$ — случайное (неизвестное) количество лет, которое вы проживете, если сейчас вам $x$ лет, $\Pr\bigl[T_x\ge t\bigr]$ — вероятность того, что $T_x$ окажется больше (желаемого) количества лет $t$ , $m$ — модальное (наиболее вероятное) значение распределения (в годах) и $b$ — дисперсия (в годах). $m$ и $b$ автор называет Гомперц-параметрами. В коде на R это выглядит так:

#' Gompertz survival probability
#' Note that the modal value isn’t the mean (average) value,
#' and that the dispersion coefficient isn’t quite the
#' standard deviation either, although both are close.
#' @param x current (conditioning) age
#' @param t time (to survive)
#' @param m modal age (mode), the age at which a person is most likely to die
#' @param b dispersion coefficient (close to SD, but not SD)
TPXG <- function(x,t,m,b) {
  exp(exp((x-m)/b)*(1-exp(t/b)))
}

Функция возвращает вероятность того, что, начиная с возраста $x$ , вам удастся протянуть ещё $t$ лет. И текущий возраст тут важен, потому что чем старше человек, тем больше его вероятность дожить до большего возраста. Повторюсь, что эта модель учитывает только старение организма, исключая всякие прочие возможности преждевременного отъезда в мир иной (для учета этих вероятностей в модели и нужно дополнение от Мейкема, однако кардинально на результат это не влияет, если, конечно, не задавать эту дополнительную вероятность очень высокой).

Вот несколько примеров (не обязательно устанавливать R на компьютер для их выполнения, можно воспользоваться каким-нибудь онлайн-сервисом):

# Вероятность дожить до 100 лет, когда тебе 0 лет, всего 1,1%.
# Здесь предполагается, что наиболее вероятная продолжительность
# жизни m = 85 лет, а дисперсия b = 10.
round(TPXG(0,100,85,10),digits=3) # → 0.011
# То же, но в возрасте 60 лет. Вероятность выросла лишь до 1,2%
round(TPXG(60,40,85,10),digits=3) # → 0.012
# Вероятность дожить до 100 лет, когда тебе уже 85 лет, 3,1%.
round(TPXG(85,15,85,10),digits=3) # → 0.031
# В книге автор приводит best-fitting m и b для мужчин в России по данным
# HMD 2011, в возрасте 60 лет вероятность дожить до 100 лет всего 0,3%.
round(TPXG(60,40,74.18,13.99),digits=3)
# Дотянуть до 80 уже вполне возможно — вероятность 31,6%
round(TPXG(60,20,74.18,13.99),digits=3)

Здесь сразу можно сделать некоторые выводы для финансового планирования. Например, что если сейчас тебе 30 лет, то рассчитывать дотянуть до 90 — это довольно консервативное планирование (с российскими $m$ и $b$ вероятность получается 4,7%), но если уже дожил до 80, то ещё 10 лет протянешь с вероятностью 20,6%, что немало. Понятно, что можно сразу допускать долгожительство вне зависимости от вероятности, но по экономической науке (см. utility) рационально всё же исходить из текущего возраста и периодически план пересматривать, исходя из актуальных вероятностей.

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

Поскольку $T_x$ — это случайная величина, на неё интересно смотреть и использовать в Монте-Карло симуляциях, для чего автор дает функцию GRAN, которая возвращает N исходов в годах дожития:

#' Generates N random conditional remaining lifetimes for an individual who is
#' currently x-years-old, based on Gompertz law
#' @param N number of samples
#' @param x current (conditioning) age
#' @param m modal age (mode), the age at which a person is most likely to die
#' @param b dispersion coefficient (close to SD, but not SD)
GRAN <- function(N,x,m,b) {
  b*log(1-log(runif(N))*exp((m-x)/b))
}

Ниже её вывод для текущего возраста 60 лет и разных $m$ и $b$ . В худшем случае прожить не удастся и года, медиана и среднее близки и сильно зависят от моды, но не равны ей, а последние долгожители заканчиваются на 110-м году жизни, причем для большего значения дисперсии $b$ при меньшей моде $m$ я стабильно получаю чуть большее максимальное значение.

> summary(60+GRAN(10000,60,85,10))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  60.02   75.01   82.39   81.76   88.94  107.39
> summary(60+GRAN(10000,60,74.18,13.99))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
  60.00   68.30   74.93   75.52   82.04  109.27

Как видно, в среднем 60-летний предпенсионер проживает ещё 15–20 лет, и только верхний квартиль счастливчиков добирается до планки 22–30 лет (в зависимости от моды и дисперсии). Традиционное у финансовых консультантов значение продолжительности пенсии в 30 лет исходит из консервативной оценки, а не медианной.

Самое интересное начинается, если запускать симуляцию для двух случайных величин $L_\xi$ — продолжительность жизни портфеля при ставке изъятий $ξ$ (кси) и $T_x$ — продолжительность жизни владельца портфеля. Нас, конечно, интересует вероятность того, что $L_\xi\ge T_x$ , т. е. что деньги не закончатся раньше времени. Или обратная ей, которую называют lifetime ruin probability (LPR) и обозначают буквой $\varphi$ (фи, varphi), которую и выдает функция VARPHI.SM. Внимательные тут заметят броуновское движение со смещением для симуляции случайного блуждания и генерацию на недельном интервале. Формулы, которые у автора есть для любой из функций, и подводящие к ним не привожу, потому что мне лень писать много LaTeX. Отправляю за ними в книгу.

#' Using simulation, approximates the lifetime ruin probability when investment
#' returns are random with parameters (ν,σ), under a Gompertz (m,b) model.
#' Uses weekly returns for simulation.
VARPHI.SM <- function(N,x,m,b,xi,nu,sigma) {
  V<-c()
  wks<-round(GRAN(N,x,m,b)*52,0)+1

  for (i in 1:N) {
    sB<-sigma*sqrt(1/52)*cumsum(rnorm(wks[i]))
    Z<-exp((nu/52)*c(1:(wks[i]))+sB)
    V[i]<-sum((1/Z)/52)
  }

  sum(V>=1/xi)/length(V)
}

Я взял какой-то средний вариант для $m$ и $b$ , текущий возраст 60 лет, доходность портфеля $\nu$ = 4% (греческая буква ню) — это реальная геометрическая доходность в непрерывном времени, т. е. не арифметическая, которая была бы $\mu$ (мю), стандартное отклонение доходностей портфеля $\sigma$ = 20% и два значения для ставки изъятий — $\xi$ = 4 и $\xi$ = 3% (это constant spending по Бенгену, т. е. определяем сумму изъятия по этой ставке в начале и далее индексируем на инфляцию, которую в расчете по реальной доходности игнорируем).

round(VARPHI.SM(10000,60,80,12,0.04,0.04,0.2),digits=3) # → 0.129
round(VARPHI.SM(10000,60,80,12,0.03,0.04,0.2),digits=3) # → 0.06

С этими значениями $m$ и $b$ продолжительность жизни для 60-летнего человека в среднем выходит около 78,5 лет (из функции GRAN), что в сочетании с указанными доходностью и волатильностью дает нам вероятность преждевременного исчерпания капитала ~13% для ставки изъятий $\xi$ = 4% и ~6% для $\xi$ = 3%.

Что мне больше всего понравилось, так это то, что примерно те же значения вероятностей можно получить и вовсе без симуляций (требующих компьюта), чисто аналитически (= мгновенно). Формулы приводить не буду, я в них сам мало что понял, да и автор выход на финальную формулу опустил, отправив за подробностями в свои статьи. Там он достал уже какую-то совсем черную магию под названием moment matching и, приговаривая что-то про стохастические дифференциальные уравнения и «исчисление Ито», выдал следующий рецепт. Забавный момент: на фоне всего матана используется преобразование непрерывной доходности в арифметическую через приблизительную формулу $\mu = \nu + \sigma^2/2$ (известный лайфхак и не самый точный из возможных, кстати).

#' Computes the incomplete Gamma function
G <- function(a,c) {
  # Avoid c=0, when a=0.
  integrand <- function(t) { (t^(a-1)*exp(-t)) }
  integrate(integrand,c,Inf)$value
}

#' Computes a Gompertz random lifetime present value
a <- function(nu,x,m,b) {
  b*exp(exp((x-m)/b)+(x-m)*nu)*G(-b*nu,exp((x-m)/b))
}

#' Same as VARPHI.SM, using moment-matching techniques.
#' The moment-matching estimates for ϕ are slightly lower (by 2–3% points)
#' compared to the simulation-based results. Note that VARPHI.SM > VARPHI.MM
#' isn’t always the case, and is actually reversed at higher withdrawal rates ξ.
#' For example, when ξ = 10%, ν = 5%, and σ = 25%, the LRP under the MM
#' technique is higher by a few percentage points compared to SM.
VARPHI.MM <- function(x,m,b,xi,nu,sigma) {
  mu <- nu+(0.5)*sigma^2
  M1 <- a(mu-sigma^2,x,m,b)
  M2 <- (a(mu-sigma^2,x,m,b)-a(2*mu-3*sigma^2,x,m,b))/(mu/2-sigma^2)
  alpha <- (2*M2-M1^2)/(M2-M1^2)
  beta <- (M2-M1^2)/(M2*M1)
  # Shape is Alpha, Scale is Beta
  pgamma(xi,shape=alpha,scale=beta,lower.tail=TRUE)
}

И, должен сказать, это работает:

round(VARPHI.MM(60,80,12,0.04,0.04,0.2),digits=3) # → 0.111
round(VARPHI.MM(60,80,12,0.03,0.04,0.2),digits=3) # → 0.053

LPR действительно получается очень близкой к той, что можно извлечь из симуляции двух случайных величин, обычно ниже на 2–3%, за исключением случаев с очень высокой ставкой изъятий (пример в комментарии к функции). Ну а я теперь читаю книжку и слушаю лекции по матану, потому что после всех увиденных интегралов в полной мере осознал свою глупость.

Вы можете спросить, какой в этом всём практический смысл? Да никакой, наверное. Я читал ради того, чтобы узнать что-нибудь новое, и узнал. Заодно получил вводный курс в R. К тому же всё это просто красиво. Ещё сравнил получаемые вероятности (из функции для симуляции только $L_\xi$ со статическим возрастом, которую тут не приводил для краткости) со своим Монте-Карло симулятором при похожих параметрах (они сошлись).

Недостатком всех рецептов, которые используют ставку изъятий, я бы, конечно, назвал их ригидность. Хотя автор и делает вполне успешную попытку задать гибкие расходы в одной из функций (зависящие от результатов портфеля, но также и «рационального потребления», основанного на величине имеющейся гарантированной пенсии, longevity risk aversion и субъективной ставке дисконтирования), его подход написания кода через формальную математику в каких-то местах уже слишком ограничивает, ведь дополнительную логику изъятий, разные доходности и суммы пополнений или расходов в разные периоды времени в прилично выглядящую формулу уже не запихнешь. Так же, как и симуляцию на исторических или бутстрап-данных. В общем, на полноценный инструмент для планирования эти рецепты вряд ли потянут (эту проблему все-таки лучше решать алгоритмически), но для удовлетворения академического интереса очень хороши.