Алгоритмвычисления сверток функций распределения высокого порядка взадачах имущественного страхования и перестрахования Текст научной статьи по специальности «Математика»

МОДЕЛИРОВАНИЕ В ЭКОНОМИКЕ

УДК 368.1:519

С. Д. Голубев, Л. А. Черная, А. Г. Шухов

АЛГОРИТМ ВЫЧИСЛЕНИЯ СВЕРТОК ФУНКЦИЙ РАСПРЕДЕЛЕНИЯ ВЫСОКОГО ПОРЯДКА В ЗАДАЧАХ ИМУЩЕСТВЕННОГО СТРАХОВАНИЯ И ПЕРЕСТРАХОВАНИЯ

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

Постановка задачи. Расчет страхового тарифа в широком классе видов имущественного (non life) страхования основывается на вычислении кумулятивной функции распределения страхового ущерба, обусловленного страховыми событиями в группе застрахованных объектов и покрываемого страховой компанией, с последующим применением одного из традиционных актуарных принципов формирования премий [6]. Вычисление кумулятивной функции распределения совокупного страхового ущерба по группе однородных объектов осуществляется по известной формуле [1]

n

R(x) = po(N)h(x) + ^ pk (N)F (k)(x), (1)

k= 1

где h(x) — функция единичного скачка (функция Хевисайда) вида

!0, если x < 0,

(2)

1, если x > 0;

здесь F(k)(x) — функция распределения совокупных страховых выплат страховщика, обусловленная наступлением k страховых событий и определяемая как k-кратная свертка функции распределения F0(x) страховых выплат страховой компании (предполагается, что застрахованные объекты однородны, а страховые выплаты, обусловленные страховыми событиями в группе застрахованных объектов, являются статистически независимыми одинаково распределенными случайными величинами, каждая из которых описывается функцией распределения F0(x)); pk (N) — вероятность наступления ровно k страховых событий (k = 1, 2, ...,n) в группе принятых на страхование объек-

тов; N — число застрахованных объектов по рассматриваемому виду страхования.

Необходимое условие адекватного описания формулой (1) кумулятивной функции распределения совокупных страховых выплат состоит в том, чтобы верхний предел суммы в формуле (1) удовле-

п

творял условию рк ^) = 1 с достаточно высокой точностью. По-

к=0

следнее означает необходимость вычисления сверток функции ^0(к)(х), к = 1, 2,... ,п, что составляет основную трудность при вычислении кумулятивной функции распределения совокупных страховых выплат (1). Практический опыт применения модели аккумуляции (1) в актуарных расчетах показывает, что число п в этой формуле достигает величины 200 ... 400, то есть возникает практическая необходимость вычисления сверток функции ^о(х), зачастую заданных в графическом или табличном виде, кратности 200 . . . 400.

Кроме того, часто функция ^0(х)имеет один или несколько разрывов.

Например, если при страховании имущественных рисков используется механизм франшизы, то с ненулевой вероятностью имеет место ситуация, когда при наступлении страхового случая фактический размер страховых выплат оказывается равным нулю, то есть в точке х = 0 функция распределения вероятности страховых выплат от единичного страхового случая имеет скачок С1. Второй пример: если размер страховых выплат по договору страхования лимитирован некоторой константой С2, то функция распределения страховых выплат ^0(х), обусловленных единичным страховым событием, в точке х = С2 также претерпевает разрыв (см. рис. 1).

Итак, численная реализация модели аккумуляции (1) для вычисления функции распределения совокупных страховых выплат и определения страховой премии связана с решением следующей математической задачи: определить функцию распределения случайной величины

I

у = Е Хк, (3)

к=1

где Хк — статистически независимые одинаково распределенные случайные величины, I = 1, 2,..., п, а п определяется исходя из условия (3), причем

р {Хк <х} = ад, у к. (4)

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

w

o.s

о.б -

0 500 1000 1500 2000 *

Рис. 1. Функция распределения страхового ущерба от единичного страхового случая с учетом франшизы и лимита

Функция распределения случайной величины У может быть формально вычислена с помощью рекуррентной процедуры вида

^(1) (ж) = ^с(ж),

х

Г(к)(ж) = ! ^(к - 1)(ж - у)^о(у), к = 2, 3 ,...и. о

При к = и соответствующая свертка дает искомую функцию распределения Р(ж) случайной величины У:

р(ж) = Р {У < ж} = ^(п)(ж), (6)

где ^(п)(ж) — и-кратная свертка функции распределения (4).

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

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

0.5

1

0

0

500

1000

1500

2000

х

Рис. 2. Графики функций распределения:

1 — непрерывная; 2 — разрывная; 3 — нормированная к единице непрерывная функция

что воспользоваться предельными теоремами "в лоб" практически не удается (см. рис. 3). Тем не менее, результатами работы [2] можно воспользоваться и весьма эффективно, если удается исходную функцию F0(x) представить в виде суммы двух составляющих — непрерывной Fн(x) и ступенчатой Fс(x) (рис.2). Именно предположение о возможности представления свертываемой функции Fo(x) в виде суммы непрерывной и ступенчатой составляющих лежит в основе предлагаемого в данной работе алгоритма исчисления сверток F0n)(x) высокой кратности.

Вывод расчетных соотношений. Будем считать известным относительно функции распределения F0(x), что в точках xV, V = 1, 2,..., т, эта функция имеет разрывы величиной су, > 0, т.е.

а в остальных точках x € (0, ж) F0(x) непрерывна. (Реально число разрывов т невелико, обычно т < 2... 3). На основании выражений (7) можно построить ступенчатую составляющую функции F0(x), определив ее формулой

V=1

где Л,(-) — функция единичного скачка (2).

Непрерывную составляющую функции F0(x) определим как разность

с^ = lim[Fo(xo + е) - Fq(Xv - е)], v = 1, 2,...,m, (7)

m

(8)

Fh(x) = Fo(x) - Fc(x).

(9)

В силу свойств исходной функции распределения, Fн(x) является непрерывной неубывающей функцией. Для Fc(x) и Fн(x) имеют место следующие предельные соотношения:

lim Fc(x) = 7 Cv,

x—^oo ' J

v=1

lim Fh(x) = ан = 1 - У^ Cv,

x

v=1

m

(10)

так как lim F0(x) = 1.

x—^^o

Чтобы получить необходимые расчетные соотношения для вычисления n-кратной свертки, воспользуемся преобразованием Лапласа исходной функции распределения F0(x) и ее составляющих:

Fo(p) = e-pxFo(x)dx,

Fc(p) = e-pxFc(x)dx,

Fn(p) = e-pxFH (x)dx.

(11)

Далее воспользуемся известным свойством преобразования Лапласа [4, гл.2, п. 9], по которому свертке оригиналов в пространстве изображений соответствует произведение изображений свертываемых функций. С привлечением данного свойства рекуррентная процедура (5) в терминах изображений записывается следующим образом:

Fo(1)(p) = Fo(p),

= Fo(k-1)(p)pFo(p), к = 2,3,..., п.

Полученные рекуррентные соотношения могут быть разрешены в замкнутой форме

(12)

F0(k)(p) = pk-1 [Fo(p)

в частности, при k = n имеем

k = 1, 2,

FOn)(p)= pn-1 Fo(p)

(13)

(13a)

С другой стороны, в силу формулы (9) изображение функции распределения F0(x) может быть записано как

Fo(p) = Fн(p) + JPс(p). (14)

k

n

Подставляя выражение (14) в (13 а) и привлекая соотношение (6), запишем формулу для изображения по Лапласу искомой функции распределения F(x):

L{F(z)} = #(p) — pn—1 F0(p)

= p

n — 1

Fu (p) + Fa (p)

(15)

Выражение (15) преобразуем, воспользовавшись формулой бинома Ньютона

- *—^ Г ~ 1 k г ~

F^(p)= Pn-^ С Fu(p) Fu(p)

n—k

(16)

k=0

П!

ГДе СП = k! • (n - k)!

— биномиальные коэффициенты (принимается,

как обычно, 0! = 1).

Внесем множитель рп-1 под знак суммы и преобразуем к-е слагаемое в формуле бинома Ньютона (16) следующим образом:

ckpn-1 ад ад

nk

— Сkpk-1

- Cnp

k г _ -| n—k

^(p) pF(p) , k — 0,1, 2,..., n. (17)

Рассмотрим алгоритм вычисления оригиналов, соответствующих изображениям (17). Сначала найдем оригиналы, соответствующие предельным значениям к = 0 и к = п. При к = 0 имеем

СП0) - р-^с(р) = СП0) - рад • р-Рс(р) •... • р^ (р). (17а)

р Ь -I р V-^-'

п сомножителей

В качестве оригинала изображения ~ • р^^с(р) = ^^с(р) естественно

р

принять Fс(x).

Каждому умножению на pFс(p) в пространстве оригиналов соответствует операция свертывания вида

Zo,i+1(x) — J Zo,i(x - y)dFa (y), i — 0,1, 2,...,n - 1), 0

Zo,o(x) — Fc(x).

(18)

Предполагая, что ^г^), г = 1, 2,... ,п — 1, представляют собой или функции распределения или их составляющие, будем считать, что ^г^) = 0, если x < 0.

n

n

x

Далее воспользуемся известным свойством интеграла Стилтьеса [5, гл. 8, п. 6, стр. 253]

X

! С(ж - - а) = С(ж - а), (19)

о

где С(ж) — произвольная непрерывная функция; — функция Хе-висайда вида (2); а — точка внутри интервала интегрирования (0, ж).

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

Полагая, что такая аппроксимация выполнена для функции ^0^(ж), найдем выражение ^0^+1(ж) с привлечением формул (18) и (19). Имеем

т

^о,г+1(ж) = У^ Су£о,г(ж - жи),г = 0,1, 2,... ,п - 1,

> (20)

^0,0 (ж) = ^с(ж),

где ж^ — координаты скачков, а с^ — величины скачков ступенчатой функции ^с(ж) (8).

Циклически выполняя вычисления по формуле (20), находим функцию £0,п(ж). Возвращаясь к выражению (17) и полагая к = п, получаем

СУ-1 (р)1 п = СП • • ^н(р) • • ... • . (21) 1 -1 4-*-'

п—1 сомножителей

Нетрудно видеть, что операции умножения на р.Рн(р) в формуле (21) соответствует операция свертывания вида

X

И+1,0(ж) = У КДж - у)^(у), г = 0,1, 2,...,п - 1, (22)

0

причем Р0,0(ж) = ^н(ж). Таким образом, имеем

Кг,0(ж) = ^н(п)(ж), (23)

где ^н(п)(ж) обозначена п-кратная свертка непрерывной составляющей ^н(ж) исходной функции распределения ^0(ж).

Наконец, обратимся к вычислению оригинала к-го слагаемого в разложении (16), 0 < к < п. В силу формулы (17а) к-е слагаемое

разложения (16) имеет вид

1 k

fc^fc-i nJ

Cn Р

FH(p) • pFC(p) • pFc(p) • ... • pFC(p)

п—к сомножителей

Приведенные выше соображения относительно п-кратной свертки

к

^н (ж) позволяют утверждать, что оригинал изображения рк—1 ^н(р) есть не что иное как к-кратная свертка непрерывной составляющей ^н(ж), т.е.

-1 J Js- 1

F (p) = *?>(х). (24)

I р ^ ^ ;

Далее, последовательное умножение изображения рк—1 ^н(р) на

р^с(р) соответствует в пространстве оригиналов операции свертывания вида

X

Жм+1(ж) = у (ж - у)^(у), г = 0,1, 2,...,п - 1), (25) о

причем Жк,0(ж) = ^н(к)(ж).

Применяя формулы (20) в предположении ЖкДж) = 0 при ж < 0 Уг, получаем рекуррентные соотношения вида

т

^к,г+1(ж) = V ^Жк,г(ж - XV), г = 0,1, 2,..., п - к - 1,

¿=1 } (26)

Жк,о(ж) = ^к)(ж).

Отметим, что формулы (26) должны использоваться для к, удовлетворяющих условию 0 < к < п, или, что то же самое, когда 1 < к < п - 1. Соответствующие вычисления для предельных значений к, т.е. к = 0 и к = п, приведены выше.

Итак, оригиналом к-го слагаемого в разложении (16) служит функция ^к,п—к (ж), рассчитываемая в зависимости от значения индекса к:

— с помощью процедуры (20), если к = 0;

— по формуле (23), если к = п;

— с помощью процедуры (26), если 1 < к < п - 1. Объединяя полученные результаты, получаем следующее выражение искомой п-кратной свертки исходной функции ^0(ж):

п—1

^о(п)(ж) = фо,п(ж) + ^ СПЖк,п(ж) + сп%,0(ж). (27)

к=1

В описанной процедуре наиболее сложным элементом является вычисление к-кратной свертки непрерывной составляющей ^н(к)(ж).

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

^нс(х) = - ВД, (28)

ан

где ан определяется из предельных соотношений (10).

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

Один из фундаментальных результатов теории вероятностей [3, гл. 9, п. 48] состоит в том, что сумма достаточно большого числа независимых случайных величин при выполнении так называемого условия Линдеберга описывается в классе безгранично-делимых законов распределения. С практической точки зрения в качестве такого безгранично-делимого закона удобно выбирать нормальный закон распределения или гамма-распределение.

В случае нормального закона распределения математическое ожидание и дисперсия суммы к независимых случайных величин, каждая из которых описывается функцией распределения ^н0(х), равны

тк = кт0,

Dk = kDo, где

00 Xmin Xmin

mo = j ж^^но(ж) = j ж^^но(ж) = Xmin — j F^^dr, o o o

Xmin Xmin

OO Г Г

D0= f x2dFH0(x)-= x2dFH0(x)=Xmin-2 xFH0(x)dx — m0; o o o

(29)

>

(30)

Хт1П — минимальный корень уравнения вида ^н0(х) = 1.

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

(0, если х < 0,

—-ххак 1е рх, если х > 0, Г(ак)

сю

где Г(а) = J ха-1е-ж^х — гамма-функция.

о

Параметры распределения (31) ак и в можно найти, приравнивая выражения математического ожидания и дисперсии их значениям, определяемым формулами (29), (30), т.е.

О ,

т = кт°,

"к 1П в^ =

откуда

2

7 т0

О = к—, в=£

Функция распределения, соответствующая плотности распределения (31), записывается в виде

0, если х < 0, Я(*>(х,а*,в) = ^ Г _ ________ п (34)

(32)

(33)

Y

О

pY^y,^,ß)dy, если x> 0.

Вычисление функций (31) и (34) на компьютерах с использованием пакетов типа Майсад не представляет трудностей, поскольку выполняется с использованием встроенных функций.

Предложенный метод вычисления сверток высокой кратности (п = 200... 300) свелся к вычислению некоторой биномиальной суммы вида (27), где СП — биномиальные коэффициенты.

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

1. Во всех расчетах сверток функции -н(х) и -с(х) заменяются на их нормированные аналоги

1 F (г)

ВД ^ (ж) = -Fн(ж) и Fc(x) ^ FCo(ж) = т^, (35)

^н 1 ^н

где ан = 1 — ^^ cv; cv — величина скачка ступенчатой функции Fc(x)

m

;

v=1

в точке ж = .

2. Ключевая формула (16) записывается следующим образом:

n

F(n)(p) = pn-1 V СПak(1 — ан)п-кРко(р)Р(П-к)(р). (36)

k=0

3. Свертки, соответствующие изображениям рп-1 —о (р)?СТ-к)(р), рассчитываются по формулам, приведенным выше (отличительная особенность предлагаемой модификации в том, что все свертки (х), (х) и РП0(х) при х ^ то стремятся к единице), но коэффициенты в формуле (27) теперь имеют вид

В = СП ан (1 — ан)п-к. (37)

4. Раскрыв выражение для СП, k = 0,1, 2,... ,n, в формуле (37), получим

n!

С

= 1;

k

n k=o 0!• n!

СП-1 = „ , ^ = (38)

(k — 1)!(n — k + 1)!

n!k = Ck k

= C

к!(п — к)!(п — к + 1) п п — к + 1'

Из соотношения (38) следует рекуррентная формула для расчета ск

n

Cnk = cn-1 n—^, k = 1, 2,, n. (39)

Используя соотношения (37) и (39), построим рекуррентную формулу для ВП. Имеем

ВП-1 = СП-1аН-1(1 — а„ )п-к+1. (40)

Рис. 3. Графики функций распределения:

ИСФ — исходная; 1 — 2-х, 2 — 3-х, 3 — 5-и, 4 — 10-ти, 5 — 15-и кратные свертки

Из выражений (37) и (40) получаем

вк

(1 - «и)

n—k

n — k + 1 аи

Bk-i

n

Ck-1ak-1(1 — аи)п-к+1

k

1 — «и

откуда

k n — k + 1

вп =

аи

k

1 — аи

в

k1

k = 1, 2,..., n.

ВП найдем, используя формулу (37):

ВП = C0(1 — ан )П = (1 — ан )П

(41)

(42)

(43)

Формулы (42) и (43) позволяют рассчитать все коэффициенты ВП,

к = 0,1, 2,..., п.

Естественное ограничение на расчет по выведенным формулам возникает в случае, когда значение Вк окажется больше предельного числа 10307. Этому числу соответствует предельная кратность свертки п = 420. При необходимости исчисления сверток более высокой кратности можно вновь воспользоваться предельными теоремами, не обращая внимания на разрывность свертываемых распределений.

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

О 2 104 4 104 6 104 S104

Рис. 4. Графики функций распределения:

5 — 15-и, 6 — 30-и, 7 — 50-и, 8 — 100, 9 — 150-и, 10 — 200-кратные свертки

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

СПИСОК ЛИТЕРАТУРЫ

1. Ш т р а у б Э. Математика имущественного страхования. - М., 1997. - 150 с.

2. П е т р о в В. В. Предельные теоремы для сумм независимых случайных величин. - М.: "Наука", 1987. - 317 с.

3. Г н е д е н к о Б. В. Курс теории вероятностей. - М.: Гос. изд-во физ.-мат. лит., 1961.-406 с.

4. Д е ч Г. Руководство к практическому применению преобразования Лапласа и Z-преобразования. - М.: "Наука", 1971. - 288 с.

5. Н а т а н с о н И. П. Теория функций вещественной переменной. - М.: Гос. изд-во технико-теоретической лит., 1957. - 552 с.

6. Бауэрс Н, Гербер Х., Джонс Д., Несбитт С., Х и км а н Д ж. Актуарная математика / Пер. с англ. под ред. В.К. Малиновского. -М.: "Янус-К", 2001. - 409 с.

Статья поступила в редакцию 7.02.2005 Сергей Дмитриевич Голубев родился в 1938 г., окончил в 1962 г. Московский авиационный институт им. С. Орджоникидзе. Канд. техн. наук, главный специалист Центра имущественного страхования ОАО "Росно". Автор 62 научных работ в области автоматического управления.

S.D. Golubev (b. 1938) graduated from Moscow Aviation Institute n. a. S. Ordzhonikidze in 1962. Ph. D. (Eng.), main specialist of Center of Property Insurance of open stock company "OAO "ROSNO". Author of 62 publications in the field of automatic control.

Людмила Александровна Черная родилась в 1945 г., окончила в 1969 г. Хабаровский политехнический институт. Канд. техн. наук, доцент кафедры "Теория механизмов и машин" МГТУ им. Н.Э. Баумана, главный специалист Центра имущественного страхования ОАО "Росно". Автор 69 научных работ в области технической механики и актуарных расчетов.

L.A. Chornaya (b. 1945) graduated from the Khabarovsk Polytechnic Institute in 1969. Ph. D. (Eng.), assoc. professor of "Theory of Mechanisms and Machines" department of the Bauman Moscow State Technical University, main specialist of Center of Property Insurance of open stock company "OAO "ROSNO". Author of 69 publications in the field of technical mechanics and actuarial calculations.

Алексей Георгиевич Шухов, окончил МГУ им. М.В.Ломоносова. Канд. физ.-мат. наук, актуарий Центра актуарных расчетов страховой компании РОСНО.

A.G. Shukhov graduated from the Lomonosov Moscow State University. Ph. D. (Phys.-Math.), actuarist of Center for Actuarial Calculations of open stock company "OAO "ROSNO".

В издательстве МГТУ им.Н.Э. Баумана вышла из печати книга

Cадовская Т.Г. Анализ бизнеса. В 4 ч. Ч. 4. Организационно-экономический анализ бизнеса. Монография, 2006 г., 280 с.

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