CC BY
88
МОДЕЛИРОВАНИЕ И ПРОГНОЗИРОВАНИЕ СМЕРТНОСТИ НА ОСНОВЕ ПУАССОНОВСКОГО ПРОЦЕССА ДЛЯ ПОДДЕРЖКИ ПРИНЯТИЯ РЕШЕНИЙ ПРИ ПЛАНИРОВАНИИ БЮДЖЕТА СОЦИАЛЬНЫХ ВЫПЛАТ РЕГИОНА РОССИЙСКОЙ ФЕДЕРАЦИИ
УДК 33.338
Ольга Ильинична Джаксумбаева,
аспирант кафедры Информационных систем в экономике, Экономический факультет, Санкт-Петербургский государственный университет Тел.: (921) 389-02-45 Эл. почта: olgadzh@rambler.ru
В статье рассмотрены существенные модельные характеристики процесса смертности населения региона Российской Федерации, произведен анализ, моделирование и прогнозирование значения показателя интенсивности смертности на год вперед. Предложена схема имитационного моделирования накопленной интенсивности смертности в течение прогнозного года по месяцам на основе неоднородного процесса Пуассона. Представлен пример расчета. Также предложен метод калибровки входных данных с целью лучшего соответствия распределению Пуассона.
Ключевые слова: социально-экономическое прогнозирование, интенсивность смертности, неоднородный процесс Пуассона, имитационное моделирование случайных процессов.
Olga I. Jaksumbaeva,
Post-graduate student, the Department of Information Systems in Economy, Economic institute, Saint-Petersburg State University, Tel.: (921) 389-02-45, E-mail: olgadzh@rambler.ru
MODELLING AND PROGNOSTICS OF MORTALITY BASED ON THE POISON PROCESS FOR DECISION SUPPORT IN PLANNING A BUDGET OF SOCIAL BENEFITS FOR A REGION OF THE RUSSIAN FEDERATION
The article introduces an approach to modeling of mortality process as a crucial factor of social benefits budget planning for Russian Federation region. Crucial properties of the process are reviewed. Analysis, modeling and one year ahead planning are presented. The non-homogeneous Poisson process based algorithm for month cumulative mortality intensity simulation within a year is proposed. Sample calculation is provided. Method of input data calibration for the best fit to Poisson distribution is developed.
Keywords: socio-economic forecasting, mortality intensity, non-homogeneous Poisson process, random process simulation.
1. Введение
После принятия в 2004 году закона о монетизации льгот 122-ФЗ обязательства по осуществлению социальных выплат населению перешли в ведение регионов Российской федерации. В то же время не были нормативно закреплены методики социально-экономического прогнозирования. В результате в каждом регионе была разработана своя методика, часто основанная исключительно на денежном выражении показателей социальных выплат. Это является причиной несопоставимости бюджетных требований регионов при составлении бюджета на федеральном уровне, а также может стать причиной недооценки необходимого объема средств и, как следствие, - задержки оказания социальной помощи населению, объективно в ней нуждающемуся.
Для большинства социальных выплат именно смертность является ключевым фактором, определяющим будущее число их получателей. Как правило, социальные выплаты осуществляются ежемесячно. В случае смерти льготодержателя выплата прекращается. Таким образом, особое значение имеет распределение числа смертей получателей льгот по месяцам.
В статье предложен подход к прогнозированию смертности населения на территории региона РФ на основе пуассоновского процесса. В отечественной литературе эта идея впервые рассмотрена в работах Староверова О. В.[1].
Решен ряд задач, соответствующих разделам основного текста. В первой части представлены результаты анализа показателя накопленной за год интенсивности смертности населения региона как временного ряда с учетом ограничений малого объема выборки демографической статистики. Во второй части излагается схема моделирования смертности по месяцам в течение года на основе неоднородного пуассоновского процесса. В третьей части приведен пример моделирования на основе данной схемы. В четвертой части приведена характеристика модели пуассоновского процесса.
В качестве объекта исследования выбран показатель интенсивности смертности населения на территории региона РФ. В актуарной математике [2] используются следующее определения. Интенсивность смертности ц(х) = р(х)/(1 - F(x)), где х > 0 - возраст, р(х) - плотность распределения случайной величины продолжительности жизниX, функция распределения которой F(x) = Р(Х< х) определяет вероятность того, что человек умрет в течение х лет с момента рождения.
Выбор именно этого показателя обоснован следующими его свойствами [3]. Во-первых, для большой выборки людей (в демографической статистике это 100 000), родившихся в определенном году, интенсивность смертности определяется в первую очередь возрастом. Функционально эта зависимость выражается законом Гомперца-Мейкема:
ц(х) = A +fx), A > 0, x > 0, fx) = R*eax, R > 0, a >0
(1)
где ц(х) - интенсивность смертности в возрасте х, А - не зависящая от возраста компонента смертности (фоновая смертность, описывает смертность, не связанную с физическим износом организма, объясняет отклонения наблюдений от экспоненты), _Дх) - зависящая от возраста компонента смертности, описывается экспонентой, Я и а - неотрицательные параметры, определяющие конфигурацию этой зависимости.
Второе свойство следует из первого: средний уровень смертности на небольших промежутках времени изменяется плавно (это верно при отсутствии демографических катастроф уровня войн и пандемий). В какой-то момент с улучшением уровня жизни интенсивность смертности стабилизируется за счет исчерпания компоненты фоновой смертности.
В-третьих, существуют хорошо исследованные аналоги изучаемого явления из других областей знания. Например, интенсивность выхода из строя станков в теории надежности или интенсивность поступления заявок на телефонную станцию в теории массового обслуживания.
Все это позволяет абстрагироваться от специфики измерения показателя, сконцентрироваться на изучении данных с использованием математического аппарата временных рядов, применять универсальные методы из смежных областей знания и является основной теоретической предпосылкой данного исследования.
2. Анализ и моделирование ежегодной интенсивности смертей
В качестве эталонного региона для расчетов использована Вологодская область, т.к. для нее имеется необходимая демографическая статистика без пропусков, и отсутствуют выбросы, связанные с демографическими катастрофами. Поскольку процесс смертности носит по большей части биологический характер, различия между регионами будут меньше, чем при анализе рождаемости, на которую в значительной мере влияет социально-экономическое и информационное окружение.
Для анализа временных рядов в настоящее время общепринято использование подхода Бокса-Дженкинса и построение моделей процессов авторегрессии и скользящего среднего (далее -АРПСС, в англоязычной литературе -АММА).
Тот факт, что чисто недетерминированный, стационарный в широком смысле случайный процесс X, может быть представлен в виде:
X, = £а,Х,
j=1
7 £,-j
(3)
Xt =
t = 1, 2, 3,
(2)
т=0
При прогнозировании при помощи модели АРПСС есть два источника ошибок. Первый - это недетерминированная составляющая модели е, (белый шум). Вторая - неточность оценки коэффициентов модели, так как оценивание производится по небольшой выборке одной из множества потенциальных реализаций. Также нужно иметь в виду ошибку спецификации. В нашем случае выбор спецификации обоснован свойствами процесса смертности.
Построение линий экспоненциального сглаживания Б, = аХ, + (1 - а)Б 1, где 0 < а < 1 - сглаживающий коэффициент, освобождает нас от дополнительных предположений о спецификации модели. Хотя все же существует определенный класс моделей, для которых использование экспоненциального сглаживания является оптимальным. Это модельЛЯ1МЛ (0,1,1). В качестве компромисса используем их совместно: АРПСС для моделирования и экспоненциальное сглаживание для прогнозирования и оценки рисков при помощи построения доверительных интервалов.
Исходными данными для анализа является показатель «Число умерших в расчете на 1000 населения за год». В паспорте показателя [5] дано следующее определение: «Общий коэффициент смертности. Рассчитывается как отношение числа умерших за период, для которого вычисляется показатель, к средней численности населения за соответствующий период и умножается на коэффициент приведения результатов расчета к годовому выражению». График значений коэффициента имеет следующий вид:
Судя по графику, можно предположить нестационарность ряда. Однако имитационное моделирование реализаций процесса авторегрессии первого порядка X, = 0,9 х X,-1 + е,, , = 1, 2, 3, ..., 21, е, ~ N(0,1) со случайно выбранными начальными значениями показало, что при таких условиях получение нестационарных на вид траекторий - довольно частый результат. Этот факт объясняется коротким рядом наблюдений. Тест Дики-Фуллера при уровне значимости 10% позволяет не отвергать гипотезу о стационарности интенсивности смертей на территории региона.
Из трех полученных в результате применения подхода Бокса-Дженкинса альтернативных спецификаций с ограничением компонент модели ЛЯ1МА в первом случае на ЛЯ: X, = 0,89 х X,-1 + е,, во втором на I и МА: ДXt = е,+0,31 х е(-1, в третьем - на ЛЯ и МЛ: X, = 0,83 х х X,-1 + 0,35 х е,ч + е„ с учетом предположения о стационарности процесса смертности была выбрана последняя. Значение информационного критерия Акаике для нее чуть ниже, чем для интегрированного варианта, и чуть выше, чем для варианта без скользящего среднего (информационные критерии применяются для сравнения качества моделей с разным числом параметров на основе суммы штрафов за число параметров и сумму остатков; чем меньше значение информационного критерия, тем лучше).
Рассмотрим теперь кривую экспоненциального сглаживания (далее - ЭС). Для выбора коэффициента сглаживания использован метод поиска по сетке. Выбор основывался на критерии наилучшего прогноза последнего наблюдения и одновременно хорошего прилегания к данным. Выбрано значение сглаживающего коэффициента а,1 = 0,3, для которого сумма квадратов
где ц - математическое ожидание процесса, г,_ т - белый шум с конечными математическим ожиданием и дисперсией, _ весовой неслучайный коэффициент, был доказан Вольдом в 1938 году в его работе «Study in the analysis of Stationary Time Series». Для удобства последующих формулировок обозначим X t : = X t _ ц.
Для комбинированного процесса авторегрессии-скользящего среднего ARMA(p, q) используем следующие обозначения [4]:
Рис. 1. Общий коэффициент смертности в Вологодской области по годам
£
t-т
№6, 2012
110
— - ■ -Эксп. стлаж. 0,3
— - - - ДИ ЭС 0,3 ниж.
— - - - Прогноз ЭС 0,3
— - ■ - ДИ ЭС 0,3 верх.
— — Эксп. сглаж. 0,5 --ДИ ЭС 0,5 ниж.
— — Прогноз ЭС 0,5
— — ДИ ЭС 0,5 верх.
— - —Эксп. сглаживание 0,15
— - -ДИ ЭС 0,15 ниж.
— - -Прогноз ЭС 0,15
— - -ДИ ЭС 0,15 верх.
АРПСС
ДИ АРПСС ниж. Прогноз АРПСС ДИ АРПСС верх. О Общий к-т смертности
Время, год
Рис. 2. Прогнозы общего коэффициента смертности
остатков минимальна с учетом точности предсказания последнего наблюдения. Дополнительно используем прилегающие к нему сверху и снизу варианты:
• а2 = а1 х 50% = 0,15 (более сильное сглаживание),
• а3 = а1 х 150% = 0,45 (более компромиссный вариант между учетом последнего наблюдения и вычисленного значения функции экспоненциального сглаживания).
На графике представлены вышеописанные модели, их прогнозы на 3 года вперед и доверительные интервалы 95%. Видно, что кривые экспоненциального сглаживания дают достаточно близкие друг к другу прогнозы и доверительные интервалы (ДИ). Доверительный интервал модели АРПСС сначала значительно уже, затем достигает границ доверительных интервалов экспоненциального сглаживания.
После того, как составлен прогноз накопленной интенсивности смертности за год, приступаем к следующему этапу - моделированию числа смертей по месяцам на основе пуассоновского процесса.
3. Схема моделирования числа смертей по месяцам на основе пуассоновского процесса
Интенсивность смертности на территории региона является пуассо-новской считающей мерой на прямом произведении двух мер: времени и меры народонаселения. На этом основываются все дальнейшие преобразо-
вания масштабов времени и народонаселения. Далее предполагается однородность по мере народонаселения.
Шаг 1. Калибровка входных данных.
- Время приводится к размерности года / е [0, 1]. Эта размерность является стандартной для финансовых вычислений, к классу которых относится задача составления бюджета социальных выплат.
- Прогноз накопленной за год интенсивность смертности на 1000 человек постоянного населения региона обозначим Л(1). С целью лучшего соответствия распределению Пуассона Л(1) приводится к другому масштабу по считающей мере народонаселения: А (1) человек за год на N =А (1) тысяч человек постоянного населения.
Масштабированная интенсивность А (1) подбирается по сетке значений на основе результатов имитационного моделирования множества реализаций пуассоновского процесса по алгоритму, приведенному в шаге 3 данной схемы. Критерий выбора основан на факте малой вероятности того, что за месяц на территории региона не умрет ни одного человека. Обозначим эту вероятность р0. Полученное значение А(1) проверяется на соответствие значению А(1/12) х КП^е^, где КП-критерий - коэффициент приведения выборочной средней А(1/12) по выборке значений интенсивности смертности по месяцам к такой размерности, для которой наилучшим образом выполняется критерий согласия Пирсона на
соответствие распределению Пуассона. Также проверяется визуальное соответствие размаха моделируемых и исторических траекторий.
Шаг 2. Выбор эталонного вектора накопленной по месяцам интенсивности смертности.
Если считать интенсивность смертности в разные месяцы постоянной, то можно использовать только данные годового накопленного коэффициента смертности А(1). Иначе предполагаем неоднородность интенсивности смертности по времени, и используем данные оперативной демографической статистики общего коэффициента смертности за последний год (А(1/12, А(2/12), ..., А(1)). Выбранный вектор также масштабируется до подобранной на шаге 1 размерности:
(А(1/12, А(2/12), ..., А(1)).
Шаг 3. Моделирование новой траектории.
- Применяем стандартный метод моделирования неоднородного пуас-соновского процесса через преобразование времени стандартизованного пуассоновского процесса с единичной интенсивностью [6].
По оси накопленной интенсивности смертности моделируются промежутки между скачками (т,), ] = 1, 2, ..., п. Соответствующие им случайные величины независимы и имеют показательное распределение F(т) = 1 - т > 0 с параметром I= 1. Моделирование происходит до тех пор, пока не будет исчерпана масштабированная прогнозная годовая интенсивность смертности т ^ + т2 + ... + тп =А (1). Предполагается, что на промежутке т умирает 1 человек на N постоянного населения региона.
- Для полученных таким образом скачков накопленной интенсивности смертности т^т + т2,т + т2 + ... тп определяем координаты по оси времени на основе ломаной линии, полученной при помощи линейной интерполяции соседних точек эталонного вектора накопленной смертности по месяцам (А(1/12; 1/12), (А(2/12); 2/12), ..., (А (1); 1). Подсчитываем число попаданий моментов смертей в каждый месяц и строим модельную траекторию накопленной интенсивности смертности.
4. Пример моделирования числа смертей по месяцам
Для расчета использован показатель «Число зарегистрированных умерших в расчете на 1000 населения (оперативные данные)».
Рис. 3. Интенсивность смертности по месяцам
Рис. 4. График сглаженной эмпирической плотности распределения ежемесячной интенсивности смертей
В паспорте показателя содержится следующее описание: «Общий коэффициент смертности. Рассчитывается как отношение числа умерших за период, для которого вычисляется показатель, к средней численности населения за соответствующий период и умножается на коэффициент приведения результатов расчета к годовому выражению» [5].
После приведения к месячному выражению его график имеет следующий вид (рис. 3).
Тест Дики-Фуллера свидетельствует о стационарности ряда. График сглаженной эмпирической плотности распределения построен в пакете Я, использован параметр сглаживания по умолчанию.
Визуально распределение смертности по месяцам соответствует нормальному закону. Гипотеза о нормальном распределении числа смертей за
месяц принимается с уровнем значимости 6% на основе критерия согласия Шапиро-Уилка на соответствие наблюдений нормальному распределению.
Обобщим приведенную выше схему моделирования следующим образом: будем моделировать число смертей в моменты скачков пуассо-новского процесса (т,),у = 1, 2, ..., п на основе нормального распределения. Параметры нормального распределения оцениваются по рассмотренной выше выборке и масштабируются для получения среднего значения равного 1. Таким образом, получаем маркированный пуассоновский процесс, в котором моменты скачков маркируются весами, равными количеству смертей, разыгранному с использованием нормального распределения. Разыгрываемые веса независимы от ведущего пуассоновского процесса и совокупно между собой.
Наилучшее согласие с пуассонов-ским распределением для значений общего коэффициента смертности по месяцам было достигнуто при Л(1/12) х КП^Ш = 136 человек на 100 000 постоянного населения в месяц. Значение параметра Л (1) было подобрано по сетке от 30 до 150 с шагом 5 на основе результатов имитационного моделирования 1000 траекторий накопленной по месяцам интенсивности смертности. Для р0 = 0,1% получено Л(1) = 95 человек на 6000 постоянного населения в год (72 000 в месяц). Различия междуЛ (1) и Л(1/12) х КПкритерий объясняются коротким рядом данных и небольшим значением р0. Остановимся на полученном значении Л (1), так как данных для более точной оценки недостаточно.
Для оценки разброса числа смертей по месяцам построено множество модельных траекторий и коридор на основе выборочного стандартного от-
120
100
I -Е
<и Б
I ?
ос а
ш Ф
£ I
щ
80
60
40
20
123456789 10 11 12 Время, мес.
-Исходная
смертность
о Смертность модель
- - - Исходная
смертность +2 ст. откл.
- - - Исходная
смертность -2 ст.откл.
Рис. 5. Результаты моделирования траекторий смертности по месяцам
0
№6, 2012
112
клонения накопленной интенсивности смертей на каждый месяц. Видно, что разброс вокруг эталонной траектории к концу года увеличивается (растет неопределенность при увеличении горизонта прогноза).
При появлении новых оперативных данных следует строить корректирующий прогноз с месяца, следующего за месяцем последнего наблюдения, и учитывающий фактически накопленную с начала года интенсивность смертности.
Оценим воздействие обычной эпидемии (не пандемии) или техногенной катастрофы на прогноз смертности. Для региона с достаточно большим населением влияние будет меньше, чем для региона с небольшим населением. Предположим, в результате падения самолета погибло 120 человек населения Вологодской области. С учетом общего населения порядка 1200000 человек и годовой интенсивности смертей порядка 16 человек на 1000 населения увеличение месячной интенсивности смертности будет незначительно - около 0,8%, и им можно пренебречь.
5. Пуассоновские процессы
Пуассоновские процессы обладают рядом полезных свойств [7], отличающих их от других точечных процессов и позволяющих использовать их в математических моделях процессов смертности населения региона.
Первое свойство - стационарность. Распределение числа событий внутри интервала времени зависит только от длины этого интервала. Вероятность того, что за интервал времени [,0, ,0 + т], т > 0 произойдет т = 0, 1, ... событий определяется как
P =
m
m!
(7)
где a = Хт в случае X = const. Параметр X > 0 _ интенсивность процесса. Интенсивность может быть как постоянной, так и переменной, то есть, зависящей от времени.
Интенсивность смертности на территории региона обусловлена половозрастной структурой населения, которая на коротких промежутках времени меняется медленно. Следовательно, количество зарегистрированных смертей зависит только от периода наблюдения и изменений
интенсивности смертности, происходящих внутри него.
Второе важное свойство - независимость приращений на непересекающихся интервалах. Количество умерших в одном месяце не влияет на количество умерших в другом месяце. То же самое относится и к непересекающимся группам наблюдаемых лиц.
Третье - ординарность. Вероятность того, что для малого интервала длины к произойдет более одного события равна о(к) - величине, которая убывает быстрее чем к. Исключение составляют несчастные случаи и события катастрофического характера, которые происходят относительно редко, и ими можно пренебречь
Промежутки времени между реализациями пуассоновского процесса ть т2, ..., тп также являются независимыми случайными величинами и распределены по показательному закону Е(т) = 1 - е т > 0, а совместное распределение моментов появления ровно п событий на фиксированном интервале [0, Т] имеет равномерное распределение.
6. Выводы
Использование пуассоновского процесса для моделирования смертности по месяцам позволяет дополнить точечный прогноз смертности (на год) вариантами реализации внутри горизонта прогноза (по месяцам), на основании которых можно осуществлять оперативное планирование бюджета социальных выплат.
Предложенная схема моделирования распределения смертей по месяцам удобна тем, что требует минимальной работы аналитика. Задача выбора масштаба интенсивности смертности решается при помощи метода Монте-Карло на основе выбранной вероятности редкого события.
Возможна доработка предложенной схемы в плане декомпозиции прогноза по полу и возрасту с учетом их особенностей.
Требуются дополнительные исследования влияния масштаба входных данных и используемых законов распределения числа смертей на результаты моделирования. Следует произвести поиск и разработку альтернативных способов построения коридора прогноза.
Литература
1. Староверов О. В. Азы математической демографии _ М.: Наука, 1997.
2. Кудрявцев А. А. Демографические основы страхования жизни. СПб: Институт страхования, 1996. _ 237 с.
3. Гаврилов Л. А., Гаврилова Н. С. Биология продолжительности жизни. _ М.: Наука, 1991.
4. Тюрин Ю. Н., Макаров А. А. Статистический анализ данных на компьютере / Под ред. В. Э. Фигурнова _ М.: ИНФРА-М, 1998. _ 528 с.
5. Демографический ежегодник России. 2010: Стат. сб./ Росстат. _ M., 2010. _ 525 с.
6. Королев В. Ю., Бенинг В. Е., Шоргин С. Я. Математические основы теории риска. _ М: ФИЗМАТЛИТ,
2007. _ 544 с.
7. Феллер В. Введение в теорию вероятностей и ее приложения, Т.1, -М.:Мир, 1984.
8. Русаков О. В. Суммы независимых пуассоновских субординаторов и их связь со строго «альфа»-устойчивы-ми процессами типа Орнштейна-Улен-бека // Записки научных семинаров ПОМИ. _ 2008. _ Т. 361. _ С. 123-127.
References
1. Staroverov O. V. Foundations of mathematical demography _ M.: Nauka,
1997.
2. Kudryavcev A. A. Demographic foundations of life inshurance. St. Petersburg: Institute of Insurance, 1996. _ 237 p.
3. Gavrilov L.A., Gavrilova N.S. The biology of life span. _ M.: Nauka, 1991.
4. Tyurin U. N., Makarov A. A. Computer statistical analysis of the data / Edit by. V. E. Figurnov _ M.: INFRA-M,
1998. _ 528 p.
5. Demographic yearbook of Russia. 2010: Stat. collection. / Rosstat. _ M., 2010. _ 525 p.
6. Korolev V. U., Bening V. E., Shorgin S. Y. Matematical basics of risk theory. _ M: PHISMATLIT, 2007. _ 544 p.
7. Feller U. An introduction to probability theory and its applications, Vol.1, -M.:Mir, 1984.
8. Rusakov O. V Sums of independent Poisson subordinators and their connections with strictly alpha-stable processes of the Ornstein-Uhlenbeck type // Notes of scientific seminars POMI. _
2008. _ Vol. 361. _ P. 123-127.
m —a
a e
