Реконструкция и исследование датчика псевдослучайных чисел в VBA-подсистеме Microsoft Office

Бородин Андрей Викторович кандидат экономических наук профессор, кафедра информатики и системного программирования, ФГБОУ ВПО «Поволжский государственный технологический университет» 424000, Россия, Республика Марий Эл, г. Йошкар-Ола, пл. Ленина, 3 Borodin Andrey Viktorovich PhD in Economics Professor, Department of Computer Science and System Programming, Volga State University of Technology 424000, Russia, respublika Marii El, g. Ioshkar-Ola, pl. Lenina, 3 bor@mari-el.com Другие публикации этого автора

Случайные числа играют огромную роль как минимум в нескольких отраслях прикладной математики. В первую очередь – это криптография. Один из канонических классов случайности, равномерная случайность, в сочетании с определенными требованиями к процедуре шифрования, а именно совпадение или преобладание мощности алфавита ключа над мощностью алфавита сообщения в условиях использования одноразовых ключей, обеспечивают существование теоретически стойких криптосистем ^[23]. Трудно переоценить роль источников случайности и в проблематике криптографических протоколов ^[24]. Другое важнейшее приложение генераторов случайных чисел в прикладной математике – это широкий класс численных методов, который очень условно можно разделить на два подкласса: «методы случайного поиска» и «методы Монте-Карло». Методы случайного поиска на современном уровне развития вычислительной техники оказываются порой единственным эффективным способом решения некоторых высокоразмерных многоэкстремальных оптимизационных задач ^[5]. Ряд задач обучения нейронных сетей также может быть описан в терминах случайного поиска ^[14]. Класс методов, объединяемых термином «методы Монте-Карло», незаменим при решении задач нахождения значений интегралов высокой и сверхвысокой кратности. Это огромный диапазон различных приложений: от теории переноса до численных экспериментов в такой абстрактной области теоретической физики, как теория струн. Финансовый инжиниринг (портфельное моделирование, сценарный анализ, стресс-тестирование) – другой утилизатор методологии Монте-Карло ^{[1, 4]}.

В то же время получение истинно случайных чисел в вычислительных устройствах с «чистой» архитектурой фон Неймана не возможно. Использование младших разрядов часов с осциллятором, независимым от тактового генератора процессора, или некоторой непредсказуемости взаимодействия человека с устройствами ввода информации (клавиатура, мышь и т. п.) в качестве источника случайности в названных выше приложениях не эффективно. Качество случайности здесь не велико, а быстродействие источников очень низкое. В связи с этим, для получения случайных чисел высокого качества используются, как правило, специальные аппаратные устройства, использующие физические случайные процессы с высокостабильными характеристиками случайности. Такие устройства дороги и стоимость их быстро возрастает с ростом быстродействия. Для преодоления этого неудобства широко используются генераторы псевдослучайных чисел (ГПСЧ). Практически все современные системы программирования содержат встроенные ГПСЧ. ГПСЧ порождают последовательности чисел, которые похожи на случайные по своим статистическим свойствам. Однако случайными эти последовательности не являются – они периодические, хотя период повторения может быть очень большим.

Каковы границы применимости того или иного ГПСЧ вместо источника истинно случайных чисел при решении той или иной задачи? Этот вопрос является вопросом номер один в отмеченных разделах прикладной математики.

Подсистема программирования Visual Basic for Application (VBA) системы автоматизации офисной деятельности Microsoft Office является основой огромного количества программных средств моделирования. Речь идет о программных продуктах в области финансового инжиниринга, в сфере оптимизации технических решений по критерию совокупной стоимости владения, в задачах конкурентного выбора дефицитного ресурса и т. п. Многие из этих продуктов используют методологию случайного поиска и методы Монте-Карло. С другой стороны, в документации на пакет Microsoft Office практически отсутствует информация о ГПСЧ, реализованном в соответствующей системе программирования. Не добавляют ясности в этот вопрос и статьи службы поддержки компании Microsoft, см., например, описание функции RAND ^[10], после прочтения которого остается больше вопросов, чем ответов. В связи с выше изложенным, в качестве цели данной работы было выбрано исследование практических аспектов использования ГПСЧ подсистемы VBA Microsoft Office в практике применения метода Монте-Карло в задачах финансового анализа. Другая цель работы – оценка допустимости применения обсуждаемого ГПСЧ для нужд элементарных криптографических приложений. Еще одной целью данной работы является реконструкция соответствующего ГПСЧ в терминах целочисленной арифметики. Последняя цель важна с позиций переноса огромного багажа теоретических знаний, накопленных в области линейных конгруэнтных методов генерации псевдослучайных чисел, на конкретную реализацию ГПСЧ.

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

Процесс управления кредитным риском включает в себя качественный и количественный аспекты ^[15]. Качественный аспект заключается в определении кредитоспособности (надежности) заемщика или контрагента. Современный подход к количественной оценке кредитного риска основывается на концепции «Value at Risk» ^{[9, 16, 17]}, ставшей общепринятым стандартом для оценки рыночных рисков. Применение данного подхода к оценке риска на уровне кредитного портфеля предполагает построение распределения случайной величины финансового результата по портфелю на каком-либо горизонте планирования и вычисления соответствующего квантиля распределения. Вычисление квантиля и есть измерение риска, ассоциированного со случайностью величины финансового результата.

В основу подхода к моделированию кредитного портфеля банка традиционно была положена идея независимого создания сценарных моделей отдельных кредитных договоров с последующим совместным анализом произвольной совокупности таких моделей ^[2].

При анализе кредитного риска большинства ссудных портфелей удобно использовать модель типового кредитного договора, включающую следующие четыре сценария ^{[11, 13]}:

1) течение кредитного договора без осложнений;

2) однократное или многократное возникновение устраненной просроченной задолженности;

3) сценарий успеха процедуры реализации залога в случае неустраненной просроченной задолженности;

4) списание безнадежной ссудной задолженности.

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

Предположим, что на момент закрытия некоторой позиции по каждому из независимых активов приведенный платежный баланс представлен случайной величиной `X_i` , `i=1, 2, ... , n` . Эти случайные величины будем считать попарно независимыми. Суммарный приведенный платежный баланс на момент закрытия позиции в таком случае составит:

`S_n=sum_(i=1)^n X_i` .

Известно ^[7], что если `X_1` , `X_2` , … , `X_n` – случайные величины с конечными математическими ожиданиями, то математическое ожидание их суммы существует и равно сумме их математических ожиданий:

`E(sum_(i=1)^n X_i)=sum_(i=1)^n E(X_i)` . (1)

Отсюда сразу следует линейность зависимости математического ожидания, как одной из мер риска, от объема портфеля.

Также известно ^[7], что если `X_1` , `X_2` , … , `X_n` – случайные величины с конечными дисперсиями и при этом они попарно независимы, то дисперсия их суммы равна сумме их дисперсий (равенство Бьенэйме):

`D(sum_(i=1)^n X_i)=sum_(i=1)^n D(X_i)` . (2)

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

Для меры риска «Value at Risk», к сожалению, не известны простые соотношения, аналогичные (1) и (2). Один из традиционных подходов к оценке этой меры риска основан на использовании метода Монте-Карло ^[1]. Для случая, когда портфель однороден, иначе говоря состоит из (примерно) одинаковых инструментов, в последние годы разработан метод точной оценки меры риска «Value at Risk» ^[12]. Используя этот результат, можно объективно исследовать поведение метода Монте-Карло.

Использование метода Монте-Карло для моделирования кредитного портфеля в указанной постановке тривиально. Каждое испытание включает в себя генерацию одного равномерно распределенного на отрезке `<sup>[0, 1]</sup>` псевдослучайного числа на каждый договор. После чего осуществляется проверка попадания этого числа в один из интервалов постоянства кумулятивной вероятности для модели договора. Финансовым результатом считается платежный баланс по сценарию договора, соответствующему правому концу интервала (разрыву кумулянты). В ходе суммирования финансовых результатов по всем договорам вычисляется итоговый платежный баланс по всему портфелю. То есть формируется финансовый результат по одному испытанию метода Монте-Карло. Учитывая равномерность исходного распределения, испытания представляют выборочное распределение финансового результата по портфелю, по которому можно оценить любую меру риска, в том числе и «Value at Risk». Чем больше испытаний, тем ближе выборочное распределение к теоретическому.

В качестве примера будем исследовать однородный портфель кредитных договоров «На неотложные нужды». Средний объем одного кредита – 90 тыс. рублей. Годовая ставка – 19%. Система выплат по договору – аннуитетная постнумерандо. Вероятности реализации сценариев составляют 0.63, 0.12, 0.11 и 0.14 соответственно. Средние затраты, связанные с устранением просроченной задолженности по одному договору, составляют 5 тыс. рублей. При списании безнадежной ссудной задолженности по неисполненному кредитному договору наряду с потерями возникают дополнительные затраты, связанные с соответствующей операционной деятельностью. Эти затраты для такого договора составляют в среднем 18 тыс. рублей. Перечисленные исходные данные будем использовать во всех численных экспериментах, связанных с методом Монте-Карло.

Рассмотрим зависимость меры риска «Value at Risk» на уровне вероятности 0.01 от количества инструментов в этом портфеле на фоне испытаний метода Монте-Карло (см. рис. 1). Точные значения меры риска «Value at Risk» получены с использованием пакета прикладных программ «МультиМИР» ^[12]. По рисунку видно, что оценка квантиля на левом хвосте распределения – задача не тривиальная. Действительно, достаточно заметить, что существенно (в разы) преобладающая часть объема конуса носителей, расположенная ниже сосредоточения точек, соответствующих испытаниям, почти пуста. Это верный признак того, что для оценки меры риска «Value at Risk» с достаточной точностью понадобится весьма значительное количество испытаний.

Рис. 1. Зависимость меры риска «Value at Risk» от количества инструментов в портфеле (красная кривая) на фоне испытаний метода Монте-Карло (синие точки). Испытания ограничены конусом носителя. Конус носителя задан коричневым лучом снизу и зеленым сверху.

Рассмотрим однородный портфель, состоящий из 256 описанных выше кредитных договоров. Выбор числа 2⁸ связан с желанием обострить проблемы, возникающие при использовании метода Монте-Карло с ГПСЧ системы программирования Microsoft Office. Почему этот выбор в указанном смысле эффективен будет ясно несколько позже. В результате первого численного эксперимента была получена зависимость оценки меры риска «Value at Risk» от количества испытаний метода Монте-Карло в диапазоне от 2¹⁰ до 2¹⁶, см. рис. 2.

Рис. 2. Зависимость оценки меры риска «Value at Risk» от количества испытаний метода Монте-Карло (синяя кривая) для модельного портфеля. Светло коричневая горизонтальная линия показывает точное значение оцениваемой меры риска. Наклонный красный отрезок прямой на второй половине графика показывает наличие очевидного тренда изменения оценки.

Если первая половина графика, представленного на рис. 2, не вызывает особых вопросов, то на второй половине графика происходит что-то странное: оценка целевой меры риска постепенно достигает точного значения, а затем примерно в том же темпе отдаляется от точного значения. Этот тренд линеен, его коэффициент детерминации `R^2=0.88` . Что происходит?

Для ответа на поставленный вопрос было решено провести второй численный эксперимент, в котором правая граница количества испытаний была отодвинута до величины 2¹⁸. Данный вычислительный эксперимент потребовал существенно больший объем ресурса процессорного времени. Эта и последующие ресурсоемкие задачи, описанные в данной статье, относительно быстро были решены за счет распараллеливания на вычислительных мощностях специализированного учебного класса экономического факультета Поволжского государственного технологического университета, состоящего из 16 компьютеров с конфигурацией Core i7–2600/16Gb/1Tb/2 Gigabit Ethernet. Результаты второго численного эксперимента приведены на рис. 3.

Рис. 3. Зависимость оценки меры риска «Value at Risk» от количества испытаний метода Монте-Карло (синяя кривая) для модельного портфеля на учетверенном интервале количеств испытаний. Светло коричневая горизонтальная линия показывает точное значение оцениваемой меры риска.

Рис. 3 проливает некоторый свет на возникшую ситуацию. Совместим отрезки `[2^16, 2^17]` , `[2^17, 2^17+2^16]` и `[2^17+2^16, 2^18]` домена зависимости оценки целевой меры риска от количества испытаний и посмотрим на соответствующие части графика этой зависимости из областей II, III и IV (рис. 4). Очевидно наличие автокорреляции в исследуемой зависимости. На основе полученных данных можно выдвинуть гипотезу: период последовательности чисел, порождаемых ГПСЧ системы программирования Microsoft Office, равен `2^8xx2^16=2^24` . Повтор периода легко объясняет неизменность оценки целевой меры риска в конце I, II, III и IV областей графика на рис. 3 (или, иначе, невозможность дальнейшего улучшения оценки в областях II, III и IV). Пока остаются не ясными причины возникновения тренда во второй части графика на рис.2.

Рис. 4. Совмещенные графики зависимости оценки меры риска «Value at Risk» от количества испытаний из областей II, III и IV рисунка 3. Для наглядности графики сдвинуты по оси ординат друг относительно друга на величину, равную 4 тыс. рублей.

На первом этапе исследования ГПСЧ попробуем определить период псевдослучайной последовательности. Одновременно с этим проверим, встречается ли в этой последовательности нулевой элемент. Для этого выполним простейший код на языке программирования VBA (см. листинг 1).

Листинг 1

В результате выполнения этого кода получаем результаты, представленные в листинге 2.

Листинг 2

Таким образом, наша гипотеза подтвердилась: период последовательности чисел, порождаемых ГПСЧ системы программирования Microsoft Office, равен `2^24` .

Теперь можно выдвинуть следующую гипотезу. Можно предположить, что для генерации псевдослучайных чисел `r_i` , `i=0, 1, ...` , с носителем `<sup>[0, 1]</sup>` в VBA Microsoft Office используется линейный конгруэнтный метод (ЛКМ) ^[6], определяемый следующим реккурентным соотношением:

`X_(i+1) = | aX_i + c |_m` , `i=0, 1, ...` , (3)

где `X_o` – начальный элемент, `X_0inZZ`, `0<=X_0<m` ; `a` – множитель, `ainZZ`, `0<=a<m` ; `c` – приращение, `cinZZ`, `0<=c<m`; `m` – модуль, `m in ZZ`, `m>0`; `ZZ` – множество целых чисел; `| * |_m` – обозначение наименьшего неотрицательного вычета по модулю `m`. Числа `X_0` , `a` , `c` и `m` однозначно определяют период последовательности псевдослучайных чисел, то есть они являются параметрами ЛКМ. Далее предположим, что используется линейная конгруэнтная последовательность максимального периода с модулем

`m = 2^24` , (4)

то есть используется ЛКМ с параметрами, для которых выполняются следующие соотношения (см. условия теоремы A ^{[6, с. 29]}, адаптированные под заданный модуль ЛКМ):

`| c |_2 = 1` , (5)

`| a |_4 = 1` . (6)

Предположим также, что для приведения последовательности `{X_i}` к последовательности `{r_i}` с носителем `<sup>[0, 1]</sup>` используется операция деления:

`r_i = X_i / m` , `i=0, 1, ...` . (7)

Как мы видели выше, легко может быть найдено нулевое псевдослучайное число `r_j = 0` , порядковым номером которого является некоторое целое неотрицательное число `j` . Этим числам соответствует элемент `X_j = 0` последовательности `{X_i}` . Теперь, используя следующее псевдослучайное число и соотношения (3) и (7), можно найти приращение:

`mr_(j+1) = X_(j+1) = | c |_m = c` . (8)

Для псевдослучайного числа `r_(j+2)` с учетом соотношений (3), (7) и (8) можно записать:

`mr_(j+2) = X_(j+2) = | ac + c |_m = | (a+1) c |_m` . (9)

Заметим, что из соотношений (4) и (5) следует взаимная простота чисел `c` и `m` , что означает существование для `c` мультипликативно обратного по модулю `m` числа `| c^-1 |_m` . Умножив левую и правую части соотношения (9) на `| c^-1 |_m` в соответствующем кольце вычетов, получим:

`| mr_(j+2) (|c^-1|_m)|_m = r_(j+2) |c^-1|_m =`

`= | |(a+1)c|_m (|c^-1|_m) |_m = | (a+1)c (|c^-1|_m) |_m = | a+1 |_m` . (10)

Поскольку `r_(j+2)!= 0` , то `| a+1 |_m = a+1` . Из этого равенства и из соотношения (10) следует, что

`a = r_(j+2) |c^-1|_m -1` . (11)

Используя для вычисления мультипликативно обратного по модулю числа обобщенный алгоритм Евклида ^[3] и применяя формулы (8) и (11) к ГПСЧ системы программирования Microsoft Office, получаем:

`c=12820163` , (12)

`a=16598013` . (13)

Нетрудно убедится, что число (12) удовлетворяет соотношению (5), а число (13) – соотношению (6). Прямые вычисления окончательно подтверждают нашу гипотезу: ГПСЧ системы программирования Microsoft Office может быть описан в терминах ЛКМ (3) с параметрами (4), (12) и (13) в сочетании с нормировочным соотношением (7).

Отметим, что ЛКМ с параметрами (4), (12) и (13) достаточно хорошо известен и всесторонне протестирован. Обычно в литературе этот датчик ассоциируют с ГПСЧ системы программирования Microsoft Visual Basic 6.0, см., например, ^{[8, 18]}. В то же время, ассоциации данного датчика с ГПСЧ системы программирования VBA в сети Internet найдено не было!

Идея использования ЛКМ для получения ключевых последовательностей в системах поточного шифрования возникла достаточно давно. Используя ЛКМ, можно сформировать последовательность, похожую на последовательность случайных равномерно распределенных чисел. Эту последовательность можно использовать для поточного шифрования в условиях приближенного соблюдения требований к теоретически стойким криптосистемам. В результате, казалось бы, можно ожидать относительно высокую стойкость такого шифра при кардинальном снижении емкости секретного канала. Ожидаемая стойкость здесь тем выше, чем больше период псевдослучайной последовательности, а емкость секретного канала снижается с величины, не меньшей суммарного объема всех шифротекстов, до объема хранения всего лишь нескольких параметров ГПСЧ. Однако, к сожалению, сегодня известен целый ряд эффективных атак на такие системы, в том числе, не зависящих от длины периода. В качестве примера можно указать атаки с известным фрагментом шифруемого сообщения ^{[21, 22]}.

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

Исследуем статистические свойства ГПСЧ системы программирования Microsoft Office. Известно ^[6], что если

`|m|_d = 0` (14)

и

`Y_i = | X_i |_d` ,