Стохастический мир

Название	Стохастический мир
Дата	27.09.2022
Размер	2.82 Mb.
Формат файла
Имя файла	stepanov1.pdf
Тип	Документы #700302
страница	14 из 20

1 ... 10 11 12 13 14 15 16 17 ... 20

x
x
s
C
x-x
s
call option
x
x
s
P
x
s
-x
put option
Если при истечении call-контракта цена актива x, лежащего в основе опциона, меньше цены исполнения x s
, то его владельцу нет смысла по- купать актив, и такое право ничего не стоит. Если же цена актива x на рынке выше, чем x s
, владелец call-опциона получает доход, равный разности рыночной цены актива и цены исполнения x ? x s
. Это и есть стоимость опциона. Диаграммы стоимости опционов одновременно явля- ются диаграммами дохода владельца опциона в момент его исполнения.

Стохастическое общество
221
•
Найдјм общую формулу для справедливой цены европейского оп- циона в произвольный момент времени. Обычно для этого используют следующую идею. Предположим, никому достоверно не известно, вырас- тет актив или нет. Это означает, что его цена x в момент исполнения, в среднем, должна быть такой же, как и сегодня x = x
0
. Однако, в силу волатильности рынка существует распределение вероятностей W (x) то- го, что x будет иметь большее или меньшее значение, чем x
0
. Будущая стоимость call-опциона равна усреднению всех возможных реализаций цены x в момент истечения контракта:
hCi =
?
Z
0
C(x)W (x)dx =
x s
Z
0 0 W (x)dx +
?
Z
x s
(x ? x s
)W (x)dx.
(8.10)
Так как диаграмма стоимости опциона при его истечении представляет собой ломаную max(x ? x s
, 0)
, то интеграл разбивается на два, первый из которых равен нулю. Таким образом, для call- и, аналогично, для put- опционов среднее значение премий в момент истечения равно:
hCi =
?
Z
x s
(x ? x s
)W (x)dx,
hP i =
x s
Z
0
(x s
? x)W (x)dx.
(8.11)
Подстановка в (
8.11
) той или иной плотности распределения вероятно- стей будущей цены x будет давать различные формулы величины премии.
Таким образом, премия опциона должна равняться величине, которая не позволяет ни покупателю, ни продавцу в среднем получить доход.
x
x
0
x
s
C(x)
W(x)
x
x
0
x
s
P(x)
W(x)
call option
put option
При увеличении страйковой цены x s
ломаная линия внутренней стои- мости call-опциона C(x) = max(x ? x s
, 0)
сдвигается вправо. Поэтому значение интеграла и, следовательно, премия падает. При уменьшении x
s
, наоборот, в зону действия колокола плотности распределения W (x)
попадает прямая x ? x s
, и премия возрастает. Для put всј наоборот.
Аналогично, по мере приближения к дате истечения неопределјнность будущей цены актива x (ширина колокола) уменьшается. Следователь- но, премия call- и put- опционов уменьшается (при неизменном x).

222
Глава 8.
•
Понятно, что, если до истечения контракта еще есть время, стои- мость опциона будет отличаться от его внутренней стоимости (I
V
=
intrinsic value). Надбавка к ней за счјт неопределенности значения цены актива в будущем называется временной стоимостью опциона: (T
V
=
time value).
Премия = Внутренняя стоимость + Временная стоимость.
Графически это соотношение представлено двумя отрезками T
V
и I
V
на кривой цены опциона в некоторый момент до даты окончания контракта:
x
x
s
C
P
x
x
s
I
V
T
V
I
V
T
V
in-the money
out-of-the money
in-the money
out-of-the money
Говорят, что call-опцион находится в деньгах (in-the money), если те- кущая цена актива выше, чем страйковая x
0
> x s
, иначе опцион вне денег (out-of-the money). Put-опцион имеет перевернутую относительно x
s диаграмму стоимости, и его владелец зарабатывает при падении цены актива, поэтому для него терминология обратная.
•
Вычтем уравнения (
8.11
) одно из другого и введјм среднюю цену актива в будущем hxi, которую на эффективном рынке обычно пред- полагают равной текущему значению hxi = x
0
. Тогда hxi + P ? C = x s
Учјт стоимости денег вносит в это соотношение некоторые изменения.
Сформируем портфель, купив акцию за x
0
, put-опцион на неј +P , и заняв короткую позицию по call-опциону ?C. Стоимость такого порт- феля x+P (x)?C(x) изменяется c изменением цены x, однако, независи- мо от еј значения, в момент истечения контракта стоимость портфеля будет в точности равна страйковой котировке x s
. Действительно:
x + P ? C =

x + 0
? (x ? x s
) = x s
x > x s
x + (x s
? x) ? 0
= x s
x < x s
Таким образом, при формировании x+P ?C портфеля мы в будущем, че- рез время ?, гарантированно получим x s
. C учјтом стоимости денег r (l
C
31
) такой актив должен сегодня стоить x s
e
?r?
. В результате получается простое равновесное соотношение:
x
0
+ P ? C = x s
e
?r?
,
(8.12)
которое называют call-put паритетом (call-put parity):

Стохастическое общество
223
•
Чувствительность премии (цены) опциона к изменению текущей ко- тировки актива x
0
и уменьшению времени до его истечения характери- зуют следующие коэффициенты:
? =
?C
?x
0
,
? =
?
2
C
?x
2 0
,
? =
?C
?t
Зная значения ? и ?, можно оценить, насколько изменится цена опциона при небольшом изменении котировки актива: x ? x
0
. Для этого необхо- димо разложить премию в ряд Тейлора:
C(x) = C
0
+ ? · (x ? x
0
) +
1 2
? · (x ? x
0
)
2
(8.13)
Подобная нелинейная зависимость C(x) позволяет cформировать из опциона C и актива x портфель с определјнными свойствам. Купим call- опцион (дающий право на одну акцию) и, взяв в долг ? = ?C/?x
0
штук акций, продадим их на рынке (займјм короткую позицию). Стоимость такого портфеля составит ?(x) = C(x)??·x. Если цена акций увеличи- вается, то короткая позиция приносит убыток, однако стоимость опциона растјт, компенсируя его.
При малых отклонениях цены акции x от начального значения x
0
пре- мия будет меняться линейно: C(x) ? C
0
+ ? · (x ? x
0
)
, а стоимость портфеля C ? ? · x будет постоянной, так как изменение стоимости оп- циона полностью компенсируется изменением короткой позиции. Дель- та единиц проданного актива на каждый опцион называется правилом дельта-хеджирования.
Рассмотрим, как изменится стоимость такого портфеля при учјте ко- эффициента ?. Если спустя время ? после формирования портфеля цена акции изменилась с x
0
до x, то доход R(x, ?) = ?(x) ? ?(x
0
)
дельта- нейтрального портфеля можно разложить в ряд Тейлора по изменению цены и времени:
R(x, ? ) =
?
2
· (x ? x
0
)
2
? ? · ?.
Чем более сильные колебания происходят на рынке (в любую сторону!),
тем больший доход мы получаем. Однако, если цена актива не изменяет- ся (x ? x
0
), то цена опциона уменьшается, и со временем этот портфель будет приносить убыток (последнее слагаемое).
Пусть на рынке ожидаются важные новости, но, приведут ли они к росту или падению, заранее не известно. Однако, если есть уверенность в неизбежности Ё
^
существенного движения цены, можно построить ?- нейтральный портфель и получить прибыль. Подобная стратегия назы- вается вероятностным арбитражом.

224
Глава 8.
8.6 Формула Блэка-Шоулза
•
Найдјм значение премии европейского опциона в рамках модели ло- гарифмического блуждания:
x n
= x n?1
exp(r n
) = x
0
exp(r
1
+ r
2
+ ... + r n
) = x
0
exp(r).
(8.14)
Если итоговая доходность через n торговых дней является случайным числом с гауссовым распределением r = µ + ??, то распределение для цены будет логнормальным (стр.
17
):
W
L
(x) =
1
x?
?
2?
exp

?
(ln(x/x
0
) ? µ)
2 2?
2

Волатильность со временем увеличивается ? = ?
0
?
?
, где ?
0
волатиль- ность единицы времени. Если ? измеряется в долях года, то ?
0
будет годовой волатильностью доходности.
Подставляя W
L
(x)
в (
8.11
) и проводя интегрирование (l H
49
) для сред- ней цены call опциона в момент его истечения, получаем:
hCi = Ї
xF
1
?
ln
Ї
x x
s

+
?
2

? x s
F
1
?
ln
Ї
x x
s

?
?
2

,
(8.15)
где Їx = hxi = x
0
e
µ+?
2
/2
среднее значение цены, а F (z) интегральное нормальное распределение (
1.12
), стр.
16
. Аналогично находится цена премии put-опциона. В этом случае необходимо просто переставить ме- стами x s
и Їx.
Для динамики актива, лежащего в основе опциона, можно заклады- вать различные параметры сноса µ (доходности) и волатильности ?. Од- нако обычно считают, что на эффективном рынке среднее будущей цены равняется текущему значению Їx = x
0
Если учитывается стоимость заимствования (временная стоимость де- нег) (l C
31
), то необходимы некоторые изменения. Пусть актив гаранти- рованно или в среднем через время ? будет стоить hxi. Тогда сегодня с учјтом процентной ставки r он должен стоить x
0
= hxi e
?r?
. Аналогич- но для цены опциона: C = hCi e
?r?
. В результате получаем известную формулу Фишера Блэка и Майрона Шоулза:
C = x
0
F
1
?
ln
x
0
e r?
x s

+
?
2

? x s
e
?r?
F
1
?
ln
x
0
e r?
x s

?
?
2

Она была выведена в 1973 г., в год открытия централизованной площад- ки по торговле опционами в Чикаго (CBOE). Рассмотрим ещј один под- ход к еј выводу, который, к тому же, будет применим и к американским опционам.

Стохастическое общество
225
•
Если цена x актива, лежащего в основе опционного контракта, ис- пытывает логарифмическое блуждание (стр.
58
):
dx = µ x dt + ? x ?W,
то премия опциона C = C(x, t) также оказывается стохастической вели- чиной. Еј изменение, в силу леммы Ито (
2.15
, стр.
55
), равно:
dC =
?C
?t
+ µ x
?C
?x
+
?
2
x
2 2
?
2
C
?x
2

dt + ? x
?C
?x
?W.
Пусть сформирован дельта-хеджированный портфель, состоящий из вы- писанного (проданного) call-опциона с премией C и купленных ? единиц базового актива (например, акции), где ? = ?C/?x. Результирующий портфель:
?(x, t) = ? · x ? C(x, t)
(8.16)
является функцией цены актива x и текущего времени t.
Будем считать, что при малых изменениях цен коэффициент ? явля- ется постоянным. Тогда изменение стоимости портфеля равно:
d? = ? · dx ? dC.
Фактически, мы каждый момент времени переформировываем портфель,
купив ? акций. После изменения цен dx, dC выбирается новая ?, и т.д.
Подставляя выражения для dx и dC, мы получим изменение стоимо- сти портфеля, которое не зависит от стохастической переменной ?W и сноса µ. Если некоторый портфель гарантированно увеличивает свою стоимость, то на эффективном рынке этот рост должен быть эквивален- тен изменению банковского депозита ?rdt с начальной суммой ?:
d? = ?
?C
?t
+
?
2
x
2 2
?
2
C
?x
2

dt = ?rdt.
Если в правую часть подставить ? из (
8.16
) и перейти ко времени, остав- шемуся до истечения опциона ? = T ? t, то получится уравнение Блэка-
Шоулза:
?C
??
+ rC =
?
2 2
x
2
?
2
C
?x
2
+ rx
?C
?x
(8.17)
Для его решения необходимо задать начальные и граничные условия.
Именно различный их выбор приводит к отличающимся результатам для опционов европейского и американского типа.

226
Глава 8.
•
Решим (
8.17
) для опционов европейского типа. Начальные условия
при ? = 0 (точнее, конечные в момент истечения) имеют вид:
C(x, 0) = max(x ? x s
, 0).
(8.18)
В уравнении (
8.17
) стоит избавиться от множителей x при производных.
Для этого необходимо перейти к переменой y = ln x. Следующая замена
C = e
?y+??
U (y, ? )
, при подходящем выборе констант ? и ?, позволяет избавиться от члена с первой производной по y, пропорционального U.
В результате получается уравнение теплопроводности (l H
50
):
?U
??
=
?
2 2
?
2
U
?y
2
Мы видели (стр.
108
), что его частным решением является гауссиана:
P (y, ? ; y
0
) =
1
?
?
2??
exp

?
(y ? y
0
)
2 2?
2
?

Общее решение линейного уравнения получается в виде суммы частных решений, соответствующих различным значениям y
0
:
U (y, ? ) =
?
Z
??
u(y
0
)P (y, ? ; y
0
)dy
0
(8.19)
Функция P (y, ?; y
0
)
имеет единственный максимум в точке y = y
0
. Его значение P (y
0
, ? ; y
0
) = 1/?
?
2??
стремится к бесконечности при ? ? 0.
Ширина колокола P (y, ?; y
0
)
при этом стремится к нулю (? - функция
Дирака, стр.
315
). Следовательно, общее решение в начальный момент времени (при ? = 0) совпадает с функцией u(x):
U (y, 0) =
?
Z
??
u(y
0
)?(y ? y
0
)dy
0
= u(y).
Поэтому u(y) имеет смысл начального значения функции U(y, ? = 0).
С учјтом проделанных замен: U(y, ?) = e
??y???
C(e y
, ? )
начальные условия (
8.18
) выглядят следующим образом:
u(y) = U (y, 0) = e
??y max(e y
? x s
, 0).
Подставляя u(y) и гауссову плотность P (y, ?; y
0
)
в общее решение (
8.19
),
мы получим формулу Блэка-Шоулза (l H
50
)

Стохастическое общество
227
•
Для опционов американского типа ситуация сложнее. Кроме началь- ных условий (
8.18
), необходимо также учитывать граничные условия.
Американский call-опцион, в отличие от европейского, не может стоить меньше своей внутренней стоимости. В противном случае его можно ку- пить за C и, сразу исполнив по страйковой цене x s
, получить гарантиро- ванный доход (x ? x s
) ? C
. Премия европейского put-опциона при боль- ших процентных ставках и x
0
< x s
, наоборот, может опускаться ниже внутренней стоимости (l C
32
).
Поэтому в случае американского опциона при решении уравнения (
8.17
)
необходимо всј время следить за тем, чтобы премия была выше внутрен- ней стоимости. Эта процедура называется задачей со свободными огра- ничениями (free boundary problem). Обычно она решается численно на дискретной решјтке (x, ?). Другой подход использование биномиаль- ной модели эволюции цены, в которой из данного состояния возможны только два перехода вверх и вниз.
•
Волатильность доходности цены актива, лежащего в основе опциона,
является ключевым параметром, который определяет его стоимость. Раз- личают историческую и подразумеваемую (implied) волатильность. Пер- вая вычисляется на основании исторических данных, а вторую рынок за- кладывает в опционные формулы, стараясь сделать некоторый прогноз будущей волатильности. Подразумеваемая волатильность может быть восстановлена из цен различных опционов. На бирже CBOE вычисляет- ся даже индекс подразумеваемой волатильности VIX.
Обычно подразумеваемая волатильность несколько возрастает при от- клонениях цены исполнения x s
от текущей цены актива x
0
, образуя па- раболоподобную зависимость ?(x s
)
, на профессиональном жаргоне име- нуемую кривой улыбки.
•
В заключение заметим, что дифференциальное уравнение Блэка-
Шоулза и одноимјнная формула для цены опциона представляет собой модель, базирующуюся на ряде допущений. Прежде всего, предполага- ется, что распределение изменений цены является логнормальным и ста- ционарным. Если статистические параметры динамики цен (в первую очередь, волатильность) не изменяются, то доходности действительно имеют близкое к нормальному распределение. Однако реальные финан- совые рынки нестационарны, волатильность изменяется со временем, и простое логарифмическое блуждание оказывается очень грубым прибли- жением к реальности. Во-вторых, непрерывное перестраивание портфе- ля для выполнения ?-хеджирования на самом деле затруднено заметной разницей между котировками bid и ask у маркет-мейкеров.

228
Глава 8.
8.7 Кривая доходности
•
Пусть B(?, t) это стоимость бескупонной облигации (векселя) в момент времени t с датой погашения t e
, т.е. через интервал ? = t e
? t
Будем считать, что еј номинальная стоимость равна единице, и, следо- вательно, B(0, t e
) = 1
. Бескупонная облигация эквивалентна депозиту с единичной стоимостью в конце. Функция B(?, t) при этом обознача- ет сумму B(?, t) < 1, которую необходимо разместить на депозите под ставку r(?, t), чтобы через время ? его величина равнялась единице.
Процентный доход облигации в момент времени t равен:
B(?, t) = e
?r(?,t) ?
=>
r(?, t) = ?
1
?
ln B(?, t).
Функция двух аргументов r(?, t) является ставкой заимствования на срок ? = t e
?t в момент времени t. По мере приближения t к дате погаше- ния t e
стоимость облигации возрастает, стремясь к своему номинальному значению, но делает это неравномерно, так как может изменяться и став- ка.
•
В фиксированный момент времени t функция r
?
= r(?, t)
, завися- щая от времени до истечения ?, называется кривой доходности (yield curve). Она определяет временную структуру процентных ставок (term structure of interest rates). Если на рынке присутствуют облигации (депо- зиты) с различными длительностями обращения ?
1
, ?
2
, ..
, то, вычисляя их эффективные доходности, мы получим, вообще говоря, различные значе- ния процентных ставок r
1
, r
2
, ..
. Их совокупность и формирует кривую доходности r
?
Кривая доходности постоянно изменяется r
?
= r
?
(t)
. Она может сдви- гаться вверх или вниз, когда все ставки изменяются на одну величину,
или определјнным образом изгибаться. Прогнозирование формы кривой доходности является исключительно важной задачей для всех участни- ков финансовых рынков.
•
Краткосрочной процентной ставкой (short term rate) r
0
(t)
назы- вают значение процентной ставки в момент времени t с истечением де- позита тут же: r
0
(t) = r(0, t)
. Естественно, мгновенных депозитов не бывает, однако, если аппроксимировать реальные данные для значений r
?
некоторой гладкой функцией, то обычно она имеет ненулевое значе- ние в точке ? = 0. Хорошим аналогом краткосрочных ставок является рынок банковских ночных заимствований для поддержания резервных требований Национального банка.

Стохастическое общество
229
•
Рассмотрим пример простой однофакторной модели описания кри- вой доходности. В ней предполагается, что динамика цен B(t, t e
)
векселя с датой истечения t e
полностью определяется динамикой краткосрочной ставки r
0
(t) = r(0, t)
. Она является единственным фактором, задающим кривую доходности. Напомню, что, если нам известна функция двух ар- гументов B(t, t e
)
, то фактически известны и форма кривой доходности r
?
(t) = r(?, t)
, где ? = t e
? t
, и еј эволюция t.
Рассмотрим портфель, состоящий из двух бескупонных облигаций с датами погашения t
1
и t
2
. Пусть отношение суммы, на которую куплен первый вексель B
1
= B(t, t
1
)
, ко второму B
2
= B(t, t
2
)
равно коэффици- енту ?. При этом второй вексель продан (куплен в короткую). В случае банка можно рассматривать выданный кредит со сроком ?
1
= t
1
? t и
полученный депозит на ?
2
= t
2
? t
. Суммарный портфель равен:
? = B
1
? ? B
2
В рамках однофакторной модели предполагается, что стоимость обли- гаций зависит от краткосрочной процентной ставки B = B(r
0
, t ? t e
)
,
которая, в свою очередь, подчиняется стохастическому процессу:
dr
0
= µ(r
0
, t) dt + ?(r
0
, t) ?W.
В этом случае стоимость портфеля также будет случайной величиной, и в силу леммы Ито его изменение равно:
d? =
?B
1
?t
+ µ(r
0
, t)
?B
1
?r
0
+
?
2
(r
0
, t)
2
?
2
B
1
?r
2 0

dt + ?(r
0
, t)
?B
1
?r
0
?W
? ?
?B
2
?t
+ µ(r
0
, t)
?B
2
?r
0
+
?
2
(r
0
, t)
2
?
2
B
2
?r
2 0

dt ? ? ?(r
0
, t)
?B
2
?r
0
?W.
Выберем долю ? таким образом, чтобы изменение портфеля не зави- село от стохастической компоненты ?W :
?B
1
?r
0
= ?
?B
2
?r
0
(8.20)
Тогда члены, пропорциональные ?W , сократятся, и динамика портфеля окажется полностью детерминированной. Если цена некоторого безрис- кового актива (в нашем случае портфеля из двух облигаций) гарантиро- ванно изменяется на d?:
d? = r
0
(t) ? dt,
то это изменение пропорционально краткосрочной процентной ставке.

230
Глава 8.
Приравняем левые части этого соотношения и уравнения, полученного по лемме Ито, подставив значение для ? (
8.20
):
?B
1
?t
+ µ
?B
1
?r
0
+
?
2 2
?
2
B
1
?r
2 0
? r
0
B
1
?B
1
?r
0
=
?B
2
?t
+ µ
?B
2
?r
0
+
?
2 2
?
2
B
2
?r
2 0
? r
0
B
2
?B
2
?r
0
Левая часть выражения зависит от t
1
, а правая от t
2
. Обе эти даты независимы, поэтому уравнение будет выполняться, если его части рав- ны некоторой функции, которая не зависит от времени истечения обли- гации. Еј принято выбирать пропорциональной функции ?, в следующем виде ?(r
0
, t) ?(r
0
, t)
. Поэтому окончательно имеем:
?B
?t
+ µ ? ? ?
? B
?r
0
+
?
2 2
?
2
B
?r
2 0
? r
0
B = 0 ,
(8.21)
где B = B(r
0
, t, t e
)
, µ = ?(r
0
, t)
, ? = ?(r
0
, t)
и ? = ?(r
0
, t)
. Это уравнение определяется двумя функциями волатильностью процентной став- ки ?(r
0
, t)
и сносом с устранјнным риском (risk adjusted drift): µ(r
0
, t) ?
?(r
0
, t) ?(r
0
, t)
. После их задания можно определить зависимость от вре- мени облигации с произвольной датой истечения t e
, а, следовательно, и кривую доходности. Еј форма будет полностью определяться одной точ- кой текущим значением краткосрочной процентной ставки r
0
Для решения дифференциального уравнения требуется задание на- чального условия. В случае с облигацией оно выбирается в следующем виде: B(r
0
, t e
, t e
) = 1
, так как в момент истечения стоимость облигации равняется единице.
Заметим, что в приведенных выше рассуждениях функция B, вообще говоря, могла быть ценой самых разнообразных финансовых инструмен- тов, поведение которых тесно связано с поведением процентной ставки.
Например, это может быть колл-опцион на краткосрочную процентную ставку с датой истечения t e
и страйковой ценой K. Для него начальные условия будут иметь следующий вид: B(r
0
, t e
, t e
) = max(r
0
(t e
) ? K, 0)
. В
случае с опционами американского типа, кроме этого, необходимо накла- дывать граничные условия.
Несмотря на теоретический характер рассмотрения динамики кри- вой доходности, мы имеем существенно феноменологическую составляю- щую в лице неизвестных функций, являющихся коэффициентами в урав- нении (
8.21
). Рассмотрим один из примеров их выбора.

Стохастическое общество
231
•
Известная модель Васичка (Vasicek, 1977) получается, если задать стохастическую динамику для блуждания краткосрочной процентной ставки в виде процесса Орнштейна-Уленбека (стр.
60
):
dr
0
= ?? · (r
0
? ?) dt + ? ?W,
где ?, ?, ? константы модели. В рамках модели предполагается, что функция ?(r
0
, t) = ?
также является некоторой константой. В результате уравнение для цены облигации:
?B
?t
+ ? ? ? r
0
? B
?r
0
+
?
2 2
?
2
B
?r
2 0
? r
0
B = 0
(8.22)
зависит от трјх параметров модели ?, ?, ? = ?? ? ?? и начального
условия B(t e
, t e
) = 1
(l C
33
).
Перейдјм к времени ? = t e
? t
, оставшемуся до истечения облигации,
и введјм процентную ставку B(r
0
, ? ) = e
?r(r
0
,? ) ?
, которая ассоциирует- ся с данной облигацией (кривую доходности). Эта кривая определяется единственным фактором r(r
0
, ? )
краткосрочной процентной ставкой. В
результате:
?r
??
? ? ? ? r
0
? r
?r
0
+
?
2
?
2
?r
?r
0

2
?
?
2 2
?
2
r
?r
2 0
+
r ? r
0
?
= 0.
(8.23)
При решении этого уравнения необходимо учитывать начальное усло- вие r(r
0
, 0) = r
0
, имеющее смысл равенства процентной ставки еј кратко- срочному значению, когда до истечения облигации времени уже не оста- лось. Прямой подстановкой можно проверить, что решением уравнения
(
8.23
) является следующее выражение:
r(r
0
, ? ) =
1
?

r
0
b(? ) +

? ? b(? )

r
?
+
?
2 4?
b
2
(? )

,
(8.24)
где введены обозначения r
?
= ?/? ? ?
2
/2?
2
и:
b(? ) =
1
?
1 ? e
?? ?
.
(8.25)
При малых ? справедливо приближјнное соотношение b(?) ? ?.
Несложно видеть, что решение (
8.24
) удовлетворяет начальному усло- вию r(r
0
, 0) = r
0
, а при ? ? ? равняется r(r
0
, ?) = r
?
. В данной модели ультрадолгосрочная процентная ставка r
?
не зависит от текущего значе- ния краткосрочной ставки r
0
и определяется только еј стохастическими параметрами ?, ?, ? и константой ?.

232
Глава 8.

Глава 9
Компьютерное моделирование
Иногда моделирование на компьютере поведения сложных стохасти- ческих систем единственный способ их исследовать. Эта глава рассчи- тана на Читателя, который любит не только формулы, но и алгоритмы.
Программирование лишь запись на очень ограниченном и формализо- ванном английском языке чјтко фиксированной последовательности дей- ствий. Эта последовательность может быть выполнена человеком, ком- пьютером или инопланетянином. Любая из программ легко переводится на человеческий язык, однако результат будет занимать больше места и,
в силу неоднозначности естественного языка, может иметь множествен- ное толкование. Глава не предназначена для обучения программирова- нию на языке C++. Для этого лучше обратиться к любому из много- численных пособий. В то же время мы по ходу дела будем приводить некоторые особенности C++, которые помогут чтению текста программ.
233

234
Глава 9.
9.1 Основы языка C++
Компьютер оперирует с целыми (int) и вещественными (oat) чис- лами. Каждое число хранится в переменной с определјнным названием.
Чтобы опечатка в еј имени не привела к недоразумению, все переменные перечисляются до их использования:
int i, j;
// объявлены две целые переменные i , j float x, y, z;
// три вещественные переменные x , y , z
Язык различает имена, написанные прописными или строчными буква- ми. Каждая законченная мысль завершается точкой с запятой, а после двух слешей // идјт произвольный текст, являющийся комментарием программиста к листингу алгоритма. С числами можно производить все- возможные арифметические операции:
x = 3.5;
y = 2;
// задајм начальные значения x и y z = x*y;
// вычисляем z , как их произведение x = -(z+y )/3;
// x теперь принимает новое значение
Все действия алгоритма осуществляются сверху вниз и слева направо,
поэтому переменная x сначала равна 3.5, затем, после вычисления z =
7 = 3.5 ? 2
, в третьей строке изменяет свој значение на -3. Кроме ариф- метических операций, существует большое число математических функ- ций, например, sin(x), натуральный логарифм log(x), и т.д. В качестве разделителя десятичных разрядов используется точка.
При работе с целыми числами необходимо помнить о некоторых осо- бенностях. Если в арифметическом выражении присутствуют только целые числа, то и результат получится целым. Иногда это может при- вести к неожиданным эффектам. Так, 7/2 равно 3 (целая часть), а не
3.5, как было бы в случае вещественных чисел.
Необходимо также помнить, что компьютер не способен проводить вы- числения с бесконечной точностью. Для целых чисел это означает, что существует максимальное значение, выше которого будет происходить переполнение. Для 32-х разрядных компьютеров целые числа по модулю меньше, чем 2'147'483'648.
Аналогичная проблема существует и для вещественных чисел. Они ограничены как по размеру, так и по точности (количеству сохраняемых разрядов числа). Для вычислений, требующих высокой точности, лучше работать не с oat, а с double вещественными переменными, дающими
двойную точность по сравнению с oat. Мы будем для обозначения вещественных чисел использовать тип Float, который может быть или обычным вещественным, или вещественным с удвоенной точностью.

Компьютерное моделирование
235
Для этого в самом начале программы необходимо вставить строку:
typedef double Float ;
// вещественный тип ,
которая во всјм коде заменит Float на double. При желании всегда можно вернуться к более быстрому типу oat, изменив только одну эту строку.
C-подобные языки обладают рядом сокращјнных обозначений, кото- рые делают их очень лаконичными. Так, очень часто при вычислении суммы необходимо прибавлять к переменной число, а результат сложе- ния снова помещать в эту же переменную. Для сокращения этой опера- ции используется следующая запись:
x +=2;
// эквивалентно x = x + 2
x++;
// эквивалентно x = x + 1
Знак ++ в названии языка обозначает переход на один уровень выше по сравнению с его прародителем, языком C. Аналогично, существуют операторы ?=, *=, /= и ??.
Операторы увеличения ++ и уменьшения ?? числа на единицу мож- но ставить как после переменной, так и перед ней. Если оператор стоит перед переменной, то она сначала изменяется, а только затем участвует в вычислениях. Если же оператор стоит после переменной, то еј измене- ние производится в последнюю очередь. Так, если x=1, то выражение y
= (x++); приведјт к y=1, а выражение y = (++x); к значению y=2.
Аналогично и для оператора ??.
•
Настоящий алгоритм начинается при появлении операции ветвления:
if(x <2)
// если x меньше 2 ,
y=5;
// то y присвоить значение 5
Сразу после if точка с запятой не ставится, так как мысль ещј не за- кончена. Двигаясь сверху вниз по алгоритму и дойдя до такой команды,
компьютер проверяет текущее значение переменной x, и только в том случае, если она меньше двух, присваивает y значение 5. Иначе пере- менная y не изменяется и программа выполняется дальше. Условный оператор может иметь и блок иначе:
if(x <2)
// если x меньше 2 ,
y=5;
// то y присвоить 5
else {
// иначе :
y=7; z=8;
// y присвоить 7 , а z 8 .
}
Если в блоках if или else выполняется не одно, а несколько действий,
то они заключаются в фигурные скобки {..}. В качестве условий можно использовать символ равенства переменных x==y и неравенства x!=y.

236
Глава 9.
•
Некоторые действия приходится делать повторно. Для этого служат операторы цикла while и for. Оператор while выполняет набор команд,
заключјнных в фигурные скобки, до тех пор, пока истинно выражение,
стоящее в круглых скобках (подобно оператору if). Так, если необходимо просуммировать целые числа от 1 до 9, мы должны написать:
int sum = 0;
// начальное значение суммы int i = 1;
// первое число while (i <10){
// пока i <10 , выполняем всј в { . . }
sum += i;
// суммируем целые числа от 1 до 9
i++;
// увеличиваем i на единицу : i=i +1
}
Первоначально целая переменная i имеет значение 1. Внутри фигурных скобок она всј время увеличивается на единицу и добавляется в пере- менную sum. Все это происходит до тех пор, пока i меньше 10. Таким образом получается результат sum=1+2+...+9=45.
Оператор for имеет три секции: инициализация, условие, при котором цикл продолжается, и операция, производящаяся в конце каждой итера- ции. Код для суммирования целых чисел, эквивалентный приведенному выше, имеет вид:
int sum = 0;
// начальное значение суммы for(int i=1; i <10; i++) // с i =1, пока i <10 , увеличивать i sum += i;
// суммируем целые числа от 1 до 9
Если в списке операторов, которые необходимо выполнять, несколько раз используется только один оператор, фигурные скобки можно не ста- вить. Отступы перед операторами, попадающими под действие команд if, while или for, делаются для повышения читаемости кода.
Внутри циклов можно использовать команды break; или continue;.
Первая прекращает выполнение цикла и выходит из блока операторов.
Вторая, наоборот, продолжает выполнение цикла, но приводит к провер- ке, не выполняя в данном такте цикла операторы, идущие ниже continue;.
Так, уже дважды реализованное выше суммирование чисел может быть выполнено также следующим образом:
int sum = 0, i=1;
// начальное значение суммы и числа while (1){
// бесконечный цикл if(i >=10) break ;
// прекращение цикла , если i >=10
sum += i;
// суммируем целые числа от 1 до 9
i++;
// увеличиваем счјтчик на 1
}
В C++ любое целое число, отличное от нуля, считается истиной, по- этому цикл while(1) никогда не останавливается. Однако, если i больше или равно 10, оператор break его прерывает.

Компьютерное моделирование
237
•
При работе с данными их необходимо где то хранить. Для этого слу- жат массивы. Они представляют собой обычные переменные, которые в математике мы обозначаем величиной с индексом a i
, b i
. При объявлении массива в квадратных скобках задајтся количество его элементов. В мас- сив можно сразу поместить набор чисел, перечислив их через запятую в фигурных скобках:
float a [10];
// массив из 10 чисел float b[3] = {1, 3, 0};
// массив из 3?х чисел : 1 , 3 , 0
a[0] = 5*b [2];
// операции с элементами массива
Массивы в языке C++ имеют нумерацию с нуля. Поэтому первый эле- мент массива это a[0], а последний a[9] (b[0] и b[2] соответственно).
•
Результаты работы программы обычно выводятся на экран или в файл. Для этого служит функция printf. Первым еј аргументом явля- ется строка, заключјнная в двойные кавычки, в которой могут распола- гаться произвольный текст и команды, определяющие, как необходимо форматировать выводимые числа. Собственно числа, которые выводят- ся, перечисляются через запятую после строки форматирования. Коман- ды форматирования всегда начинаются со знака процента %. Их количе- ство и последовательность должны совпадать с переменными, идущими после строки форматирования:
printf ("Ku -ku\n");
// вывод строки и переход ниже printf ("i=%d,j=%d", i,j); // вывод двух целых (%d ) переменных printf ("x=%g\n", x);
// вывод вещественного (%g ) числа printf ("x =%4.2 f\n", x);
// то же, с 2?мя десятичными знаками
•
Минимальная программа, которую поймјт компилятор и которая выведет на экран (в консоль) надпись Hi Stochastic World!, имеет сле- дующий вид:
# include
// подключение библиотеки вывода void main ()
// главная процедура программы
{
printf ("Hi Stochastic World !\n");
}
Основной кусок кода void main(){...} является главной процедурой, те- ло которой, находящееся в фигурных скобках, начинает выполняться при запуске программы. Воспользовавшись бесплатными компилятора- ми bcc32.exe или gcc.exe, можно откомпилировать этот текст, помещјн- ный в файл hi.cpp, в выполняемый файл hi.exe при помощи команды
bcc32.exe hi.cpp.

238
Глава 9.
9.2 Статистики
•
Если некоторые куски кода используются несколько раз, то их стоит оформить, как процедуру или функцию. Отличие между ними состоит в том, что функция возвращает значение определјнного типа, тогда как процедура нет. Перед именем функции указывается еј тип (int, oat и т.д.), а перед именем процедуры пишется слово void, что означает пу- стота, т.е. отсутствие возвращаемого значения. Процедура или функция могут иметь список аргументов, необходимых им для проведения вычис- лений. Для каждой переменной из списка аргументов указывается еј тип и могут быть заданы значения по умолчанию. В этом случае функция может быть вызвана без этих аргументов.
Запишем программу, вычисляющую среднее значение ряда данных:
# include
// подключение библиотеки вывода typedef double Float ;
// вещественный тип
Float Aver ( Float *x, int n) // с р е д н е е массива x длиной n
{
if(n <1) return 0;
// нет данных ? выходим
Float av = 0;
// начальное значение суммы for(int i=0; iav += x[i];
// суммируем данные return av/n;
// с р е д н е е значение
}
void main ()
// главная процедура программы
{
const int n = 5;
// количество данных в массиве
Float x[n] = { -2.7 , -0.2, -1.4, -1.0, -0.5};
Float av = Aver (x, n);
printf (" average =%g\n", av );
}
Количество данных в массиве помечено словом const. Это означает, что n
не переменная, а константа, и она не может изменяться по ходу ал- горитма. Для объявления статического массива необходимо указывать его размер только при помощи константы. В цикле мы складываем все элементы массива. После оператора return в функции идјт то значение,
которое мы хотим вернуть. В нашем случае это сумма значений массива,
делјнная на их количество.
Обратим внимание, что перед именем массива x в аргументах функ- ции Aver стоит звјздочка. Это указатель на имя массива, позволяющий отличить его от обычной скалярной переменной.

Компьютерное моделирование
239
•
Далее мы будем вычислять средние значения стохастических про- цессов по конечному числу n выборочных траекторий. Выясним, какой типичной ошибки при этом можно ожидать.
По n измерениям величины x найдјм арифметическое среднее:
?
x =
x
1
+ ... + x n
n
=
1
n n
X
i=1
x i
,
(9.1)
которое помечено тильдой, чтобы не смешивать его с истинным сред- ним по бесконечному числу наблюдений Їx. Будем называть ?x выбороч- ным средним. Если повторить эксперимент, то получится другой набор измерений x i
и, следовательно, другое выборочное среднее ?x. На самом деле, результат каждого измерения x i
можно рассматривать как случай- ную переменную. Выборочное среднее, являющееся суммой таких слу- чайных величин, естественно, также оказывается случайным. Таким об- разом,
выборочное среднее это случайная величина, характеризующа- яся, в свою очередь, средним, волатильностью и плотностью распределения вероятностей.
Усреднение всех возможных выборочных средних даст нам истинное
среднее:
h?
xi =
1
n n
X
i=1
hx i
i =
n n
hxi = Ї
x.
(9.2)
Интуитивно понятно, что, чем больше мы возьмјм чисел, тем ближе вы- борочное среднее окажется к истинному. Другими словами, волатиль- ность выборочного среднего должна уменьшаться с ростом n. Действи- тельно, вычислим его дисперсию:
?
2
?
x
=
?
x
2
? h?
xi
2
=
1
n
2
n
X
i,j=1
hx i
x j
i ? Ї
x
2
Все x i
являются равноправными случайными величинами. Двойная сум- ма по i, j содержит n
2
слагаемых. Из них n с одинаковыми индексами i = j
, а остальные n
2
? n
с не совпадающими равноправными индекса- ми i 6= j. Поэтому в сумме будет n слагаемых вида x
2 1
и n
2
? n вида hx
1
x
2
i
(см. также l C
15
):
?
2
?
x
=
1
n
2
n x
2 1
+ (n
2
? n) hx
1
x
2
i
? hxi
2
Если различные измерения x i
и x j
являются независимыми случайными величинами, то hx
1
x
2
i = hx
1
i hx
2
i = hxi
2
. Соответственно, x
2 1
= x
2

240
Глава 9.
В итоге получаем формулу, связывающую волатильность выборочного среднего n чисел и дисперсию ?
2
=
x
2
? hxi
2
по всей генеральной совокупности (бесконечной выборке):
?
?
x
=
?
?
n
Таким образом, волатильность выборочного среднего убывает, как ко- рень квадратный из n. На самом деле она уменьшается достаточно мед- ленно. Если мы хотим снизить ошибку измерения выборочного среднего в 10 раз, то необходимо увеличить число измерений в 100 раз!
Типичный диапазон, в который может попасть выборочное среднее как случайная величина, имеет вид:
?
x = Ї
x ±
?
?
n
В случае нормального распределения такая запись означает, что выбо- рочное среднее, полученное по n измерениям, отклоняется от истинного среднего не более, чем на ?/
?
n с вероятностью порядка 0.68.
Естественно, что реальная ? распределения обычно неизвестна. Одна- ко еј тоже можно оценить по конечной выборке ??. В этом случае сред- нее генеральной совокупности Їx можно представить в виде диапазона с t
стандартными ошибками:
Ї
x = ?
x ± t
?
?
?
n
(9.3)
При больших n, в силу предельной теоремы, отклонение от среднего обычно является нормальной случайной величиной с волатильностью
?/
?
n
. Поэтому вероятность того, что в данной серии измерений ?x откло- нится от Їx на величину, не превышающую a = t?
n
, ?
n
= ?/
?
n
, равна:
p(Ї
x ? a 6 ?
x 6 Ї
x + a) =
Ї
x+a
Z
Ї
x?a e
?
(x?Ї
x)2 2?2
n dx
?
n
?
2?
=
a/?
n
Z
?a/?
n e
?z
2
/2
?
2?
dz = ?(a/?
n
).
Задавая те или иные значения t = 1, 2, 3 (расширяя интервал), мы можем повышать вероятность попадания истинного среднего в диапазон (
9.3
).
Часто используемые значения вероятностей ?(1) = 0.683, ?(2) = 0.955
и ?(3) = 0.997.
Если ?/Їx ? 1, то относительная ошибка измерения среднего будет равна 1/
?
n
. Чтобы достичь точности 10
?3
(2-3 значащих цифр) необхо- димо поставить миллион экспериментов (n = 1000000).

Компьютерное моделирование
241
•
Для оценки дисперсии по небольшой выборке из n чисел x i
необходи- мы определенные ухищрения. Простое арифметическое усреднение при вычислении моментов соответствующих степеней приводит к неверным результатам.
Если мы возьмјм выборку из n чисел x
1
, ..., x n
и по ним вычислим выборочную дисперсию, усредняя квадраты отклонений (x i
? ?
x)
2
от вы- борочного среднего ?x, то получится некоторое значение не совпадающее с истинной дисперсией распределения D = (x ? Їx)
2
. Этот эксперимент можно повторить большое число раз. Случайные числа будут разными, и,
аналогично выборочному среднему, разными будут выборочные диспер- сии. Желательно, чтобы после усреднения по всем экспериментальным выборочным дисперсиям получился результат, совпадающий с диспер- сией реального распределения. Это произойдјт, только если мы будем вычислять выборочную дисперсию при помощи следующей формулы:
?
D =
1
n ? 1
n
X
i=1
(x i
? ?
x)
2
,
?
x =
1
n n
X
i=1
x i
(9.4)
Заметим, что в определении выборочной дисперсии в знаменателе стоит n?1
а не n, как было бы при простом арифметическом усреднении. Такая оценка называется несмещенной, так как именно она в среднем будет приводить к истинной дисперсии D. Действительно, возведя в квадрат,
перепишем выборочную дисперсию в следующем виде:
?
D =
1
n ? 1
n
X
i=1
(x i
? ?
x)
2
=
1
n ? 1
n
X
i=1
x
2
i
?
n n ? 1
?
x
2
Усредним это выражение по большому числу реализаций возможных вы- борок (по всей генеральной совокупности):
D ?
D
E
=
1
n ? 1
n
X
i=1
x
2
i
?
1
n(n ? 1)
n
X
i,j=1
hx i
x j
i =
n ? 1
n ? 1

x
2
? hxi
2

= D.
Вычисление двойной суммы проводится также, как и при выводе ?
?
x
. По- этому, если результаты наблюдений x i
можно рассматривать как незави- симые случайные числа, то усреднение выборочной дисперсии с факто- ром 1/(n ? 1) будет приводить к истинной для всей генеральной сово- купности.
Естественно, подобная поправка существенна только для малых n, ко- гда статистическая значимость результатов и так мала. Заметим также,
что несмещјнность дисперсии ?
D
, вообще говоря, не означает несмещјн- ности волатильности ?? =
p ?
D

242
Глава 9.
•
Соберјм часто использующиеся статистические функции в файле
stat.cpp:
# include
// p r i n t f
# include
// srand
# include

// s q r t , l o g typedef double Float ;// вещественный тип//??????????????????????????????????????????????????????????inline Float Sqr( Float v){ return v*v; }// квадрат числа//??????????????????????????????????????????????????????????Float Aver ( Float *x, int n) // с р е д н е е массива x длиной n{if(n <1) return 0;// нет данных ? выходимFloat av = 0;// начальное значение суммы for(int i=0; i}//??????????????????????????????????????????????????????????Float Sigma ( Float *x, int n)// волатильность массива x{if(n <2) return 0;// нет данных ? выходимFloat av = 0, di =0;// начальное значение суммы for(int i=0; i 0? sqrt (di ):0; // волатильность}//??????????????????????????????????????????????????????????Float Max( Float *x, int n)// максимум массива x длиной n{if(n <1) return 0;// нет данных ? выходимFloat max = x [0];// начальное значение for(int i=1; i x[i]) min=x[i];return min;}

Компьютерное моделирование243// Гистограмма вероятностей p для массива x длиной n// Интервал [ min . . max ] разбит на m отрезков// Если x [ i ]max он не попадает в p ,// поэтому сумма вероятностей может быть меньше единицы//void Histo ( Float *x,int n, Float *p,int m, Float min , Float max){for(int k=0; k=0 && k

244Глава 9.9.3 Случайные числа•При проведении численных моделирований нам необходимо уметь ге- нерить случайные числа. Современные компьютеры в большинстве својм являются цифровыми, а не аналоговыми устройствами. Поэтому возмож- на только генерация псевдослучайных чисел, являющихся результатом работы некоторого детерминированного алгоритма. В основе получения случайных чисел с различными распределениями лежит генератор рав- номерно распределјнных, обычно целых чисел из диапазона [0...max].Понятно, что, если необходимо вещественное равномерно распределјн- ное число в диапазоне [0...1], то целое случайное число делится на мак- симальное значение max.В C++ есть встроенная функция rand(), возвращающая целое, равно- мерно распределјнное случайное число в диапазоне от нуля до константыRAND_MAX, которая в большинстве компиляторов равна 32767.Если в программе при помощи rand() генерить случайные числа, то при еј повторном запуске получится в точности такая же последова- тельность чисел. Чтобы еј изменить, можно в начале вызвать функцию srand(seed), где целое число seed задајт номер последовательности слу- чайных чисел. Если передать в эту функцию количество секунд, прошед- ших с 1970 г., возвращаемое функцией time(0), то последовательности будут получаться каждый раз разные.Ниже в примере происходит генерация десяти случайных веществен- ных чисел из диапазона 0 6 x 6 1:# include // для функции time# include // для функций rand , srand# include // для функции p r i n t f typedef double Float ;// вещественный тип void main ( void ){srand ( time (0) % RAND_MAX );// "встряхиваем " генератор for(int i=0; i <10; i++)printf (" %12.8 f\n", Float ( rand ()) / RAND_MAX );}Если закомментировать srand, то при последовательных вызовах про- граммы ряд случайных чисел будет получаться один и тот же. Операция% обозначает остаток от деления. Например, 47%10 равно 7. С еј помо- щью мы обрезаем значение seed величиной RAND_MAX-1. Для преобразо- вания целых чисел в вещественные используется операция Float(rand()).Напомню, что 1/2 равно 0, тогда как Float(1)/2 равно 0.5.

Компьютерное моделирование245•В основе алгоритма rand() лежит формула, подобная следующей:x t+1= (16807 · x t)mod (2 31? 1).(9.5)Начиная с некоторого целого числа x0 6= 0, задаваемого при помощи функции srand(), при каждом вызове rand() происходит вычисление но- вого псевдослучайного числа на основе предыдущего. Магические кон- станты, использующиеся в этом алгоритме, являются простыми числа- ми, и качество генератора существенно от них зависит. Реализация этого алгоритма имеет вид:unsigned int gRndSeed = 1;// последнее случайное число inline void SRnd ( unsigned int seed ) // задание gRndSeed{gRndSeed = ( seed ==0)? 1: seed ;}//??????????????????????????????????????????????????????????inline unsigned int RndI ()// равномерное сл . число{return gRndSeed = (16807 L* gRndSeed ) % 2147483647 L;}От хорошего генератора мы ожидаем выполнения четырех требований:1) случайные числа равномерно распределены;2) они независимы;3) их последовательность имеет большой период повторения;4) они различны.Последние два требования не являются эквивалентными. Например пло- хой генератор может, имея формально большое RAND_MAX, выдавать толь- ко, скажем 1000 различных чисел, переставляя их самым разнообразным образом. При этом период повторения будет велик, но при построении ве- щественного генератора из всего континуума мы будем получать только1000 случайных значений. Как мы увидим ниже, эта проблема присуща как библиотечной функции rand(), так и алгоритму (9.5).Поэтому мы будем использовать генератор случайных чисел, основан- ный на эффекте переполнения 32-разрядных целых чисел. Соответству- ющая функция имеет вид:inline unsigned int RndU ()// равномерное сл . число{return gRndSeed = gRndSeed *1664525 L + 1013904223 L;}Именно она будет в дальнейшем использоваться для генерации равно- мерно распределјнных случайных чисел.

246Глава 9.•Проведјм сравнение статистических свойств трјх генераторов rand(),RndI() и RndU(). Будем сразу работать с вещественными случайными числами в диапазоне [0...1]. Так, алгоритм переполнения RndU() имеет следующую вещественную реализацию:inline Float Rnd ()// равномерное сл . число [ 0 . . 1 ]{return Float ( RndU ()) / Float (0 xffffffff );}Константа 0xffffffff является записью максимального беззнакового целого 2 32? 1в шестнадцатеричной системе исчисления. Если целые числа генерятся при помощи rand(), делить необходимо на RAND_MAX, а для RndI() на 2147483647 = 2 31? 1В качестве статистик, характеризующих равномерность случайных чисел, мы будем использовать среднее, точное значение которого равно x = 1/2, квадрат среднего x2 = 1/3, минимальное x min= 0и макси- мальное значение x min= 1. Кроме этого, будем вычислять среднеквад- ратичное отклонение вероятности попадания в один из m одинаковых интервалов, на которые разобьјм диапазон [0..1]:z m=n mmXk=1p k?1m2,где n число сгенерјнных случайных чисел, p k= n k/n, а n k-количество чисел, попавших в k-й интервал. Мы умножили эту статистику на n, так как подобное отклонение убывает, как 1/n, поэтому при любых n она будет примерно постоянной. Зная z m, можно по ?2-критерию оценить достоверность гипотезы о равномерном распределении, однако мы будем использовать z mтолько для сравнения различных генераторов.Характеристиками независимости будут служить ковариационные ко- эффициенты для моментов случайных чисел, умноженные на?n:c jk=?n x jt xk t?1 ? x jt xk t?1В случае независимых величин чисел они должны быть равны нулю.Для вычисления количества различных случайных чисел, выдаваемых генератором, воспользуемся библиотечным классом map который хранит пары ключ-значение. В нашем случае ключом будет служить случайное число, а значением количество его появлений. Класс map сравнитель- но быстро ищет и вставляет новые числа, и его размер size сообщает о количестве различных ключей (случайных чисел).

Компьютерное моделирование247# include " stat .cpp"// статистическая библиотека# include // map [ ]map lst;const int n= 1000000;// величина выборки const int nex = 1000;// число экспериментов const int m= 100;// число интервалов void main ( void ){Float a[nex], s[nex], z[nex ];Float c11[nex], c22[nex], c12[nex], c21[nex ];Float min[nex], max[nex], size [nex ];Float p[m];for(int ex =0; ex =m) k=m -1; if(k <0) k=0;p[k ]++;if(max[ex]rnd) min[ex ]= rnd;}a[ex ]/=n; s[ex ]/=n;c11[ex ]/=n; c11[ex ]=( c11[ex]-a[ex ]*a[ex ])* sqrt ( Float (n));c12[ex ]/=n; c12[ex ]=( c12[ex]-a[ex ]*s[ex ])* sqrt ( Float (n));c21[ex ]/=n; c21[ex ]=( c21[ex]-a[ex ]*s[ex ])* sqrt ( Float (n));c22[ex ]/=n; c22[ex ]=( c22[ex]-s[ex ]*s[ex ])* sqrt ( Float (n));z[ex ]=0;for(int k=0; k

248Глава 9.•После получения на основании nex экспериментов массивов со ста- тистиками при помощи функций Aver(), Sigma(), находящихся в файле stat.cpp, вычисляются среднее значение статистики и одна стандартная ошибка: Sigma(a, nex)/sqrt(nex). Ниже в таблицах стандартная ошибка используется для вывода только значащих цифр в статистике и указы- вается в виде индекса. Так, например, 0.025 3означает, что в последнем знаке может быть ошибка, т.е. 0.025 ± 0.003. Приведјм результаты рабо- ты программы:11 12 21 22size/nRndU 0.003 30.004 30.003 30.003 31.000000rand0.001 30.001 30.002 30.001 30.032767RndI 0.25 20.27 30.27 30.28 30.049Корреляционные свойства RndU и rand примерно одинаковы, а для RndI заметно хуже. Генератор RndU существенно лучше по количеству раз- личных чисел, которые он генерит. Так, на выборке из n = 1000000он выдајт все числа различными, тогда как для rand различных чисел только RAND_MAX=32767. Аналогичная проблема и у RndI, к тому же их количество существенно зависит от начального значения gRndSeed.Если нам необходимо равномерно заполнить некоторый отрезок точка- ми (например, при Монте-Карло вычислении интеграла), то, используя rand или RndI, мы получим около 10 5точек, и проведение большего чис- ла экспериментов становится бессмысленным.Статистики равномерности имеют следующие значения:x x2z100min maxRndU 0.49999 10.33332 10.00989 50.000000 0.9999999rand0.49999 10.33332 10.01009 50.000000 1.0000000RndI 0.5004 10.3337 10.73 10.000019 10.9999859 2Средние значения у первых двух генераторов в пределах стандартной ошибки совпадают с точными значениями. Генератор RndI хуже как по этим показателям, так и по степени равномерности для m = 100 интер- валов z100В дальнейшем для равномерно распределјнных чисел мы будем ис- пользовать RndU(). С помощью RndU(), а точнее, еј вещественного ана- лога Rnd(), можно получать и нормально распределјнные числа, кото- рые, собственно, нам и необходимы для стохастических дифференциаль- ных уравнений. Рассмотрим соответствующий алгоритм.

Компьютерное моделирование249•Для генерации нормально распределјнного случайного числа с нуле- вым средним и единичной волатильностью служит метод Бокса-Миллера.Перейдјм от равномерно распределјнных на интервале [0..1] чисел x1и x2к переменным: y1=??2 ln x1cos(2?x2), y2=??2 ln x1sin(2?x2). Ко- личество чисел, попавших в некоторую область двухмерного простран- ства ?, определяется двойным интегралом. При замене переменных в интеграле объјм умножается на якобиан:Z?dx1dx2=Z??x1?y1?x1?y2?x2?y1?x2?y2dy1dy2=Z?1?2?e?y2 12 1?2?e?y2 22dy1dy2Таким образом, переменные y1и y2являются независимыми (интегралы расщепляются) и имеют нормальные распределения.Функции cos и sin вычисляются достаточно долго. Чтобы избежать работы с ними, воспользуемся следующим пријмом. Выберем область интегрирования в виде круга с единичным радиусом и центром в начале координат. Будем случайным образом равномерно заполнять его точка- ми. Пусть v1и v2 координаты точки внутри круга v2 1+ v2 2= r2< 1В полярных координатах v1= r cos ?, v2= r sin ?. Так как v1и v2распределены равномерно, то и полярный угол ? будет равномерно рас- пределјн в интервале [0..2?]. Поэтому для вычисления значения cos и sin можно воспользоваться отношениями v1/r и v2/r:Float RndG ()// г а у с с о в о сл . число{static int was = 0;// была ли вычислена пара чисел static Float r= 0;// предыдущее случ . число из пары if(was ){ was =0; return r;} // r из предыдущих вычисленийFloat s, v1 ,v2;do{// ждјм попадания в круг v1 = 2* Rnd () -1;// точка в квадрате [ ? 1 . . 1 ]v2 = 2* Rnd () -1;s=v1*v1 + v2*v2;// квадрат расстояния от центра} while (s >=1.0 || s ==0);was= 1;r= v1* sqrt ( -2* log(s)/s); // первое число ( запомнили )return v2* sqrt ( -2* log(s)/s); // второе число ( вернули )}Переменные was и r объявлены в функции RndG, как статические (static).По завершении работы функции их значения будут сохранены, и при сле- дующем вызове мы будем иметь значение одного из случайных гауссовых чисел, вычисленных на предыдущей итерации.

250Глава 9.9.4 Моделирование стохастических процессовРассмотрим алгоритм моделирования одномерного стохастического про- цесса:dx = a(x, t)dt + b(x, t)?W.Если задано начальное условие x0= x(0), то последующие значения процесса x k= x(t k)получаются при помощи итерационной схемы (стр.49):x k+1= x k+ a(x k, t k) ?t + b(x k, t k)??t ?kЧем меньше временной шаг ?t, тем ближе будут свойства получаемого выборочного процесса к точному результату:# include " stat .cpp"// файл с RndG ( )inline Float a( Float x, Float t){ return 0; } // снос inline Float b( Float x, Float t){ return x; } // волат .void main (){FILE *out = fopen ("ito.out", "w");SRnd ( time (0));// "встряхиваем " генераторFloat step = 0.1;// шаг по времени для табуляции x const int num = 100; // число точек табуляции int lag = 1000;// количество дроблений шага s t e pFloat dt= step /lag , sqrt_dt = sqrt (dt );Float x=1;// начальное значение xFloat t=0;// в момент времени t fprintf (out , "%g\t%g\n", t, x);for(int k=1; k <= num; k ++){for(int j=0; j

Компьютерное моделирование251Далее идут описания функций сноса и волатильности. В примере вы- ше выбран процесс логарифмического блуждания с нулевым сносом и единичной волатильностью (стр.58).В главной функции main открывается файл с именем ito.out. В него будут помещаться результаты работы программы. Для этого объявля- ется переменная (указатель) out типа FILE, значение которой задајтся функцией fopen. Последний еј аргумент w определяет способ открытия файла. В данном случае это будет вновь создаваемый файл, открываю- щийся для записи.Функция SRnd встряхивает генератор случайных чисел. В результа- те при каждом запуске программы будет получаться новая реализация случайного процесса.В программе вычисляется num = 100 значений случайного процес- са x(t), с равным шагом по времени: t = 0, step, 2 · step, ..., num · step.Каждый интервал между вычисляемыми значениями x разбивается на большое число lag точек, время между которыми равно dt = step/lag.

1 ... 10 11 12 13 14 15 16 17 ... 20