Научная статья на тему 'Диффузионная аппроксимация математической модели деятельности некоммерческого фонда при релейном управлении капиталом'

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

CC BY
99
40
i Надоели баннеры? Вы всегда можете отключить рекламу.

Аннотация научной статьи по математике, автор научной работы — Лившиц Климентий Исаакович, Шифердекер Ирина Юрьевна

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

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по математике , автор научной работы — Лившиц Климентий Исаакович, Шифердекер Ирина Юрьевна

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Diffuse approximation of mathematical model of incomercial fund functioning under the relay control of its capital

The statistical characteristics of functioning of an incomercial fund under the diffuse approximation of its model and the relay control of the fund capital are investigated

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

К.И. Лившиц, И.Ю. Шифердекер

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

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

Математическая модель изменения капитала фонда

Как и в [1], под некоммерческим фондом понимается организация, созданная для сбора и распределения денежных средств без получения прибыли. Основной характеристикой состояния фонда является его капитал s(t) в момент времени t. В работе предполагается, что с капиталом s(t) могут происходить следующие изменения:

1. В фонд поступают денежные средства. Будем считать, что моменты поступления денежных средств образуют пуассоновский поток с интенсивностью X. Поступающие денежные суммы (премии) являются независимыми одинаково распределенными величинами с плотностью распределения ф(х), средним значением M{х} = а и вторым моментом M {х2} = а2.

2. Фонд расходует поступившие средства. Будем считать, что расходование денежных средств происходит непрерывно во времени со скоростью b(s), так что за время At выплата составляет b(s)At. Предполагается, что управление расходованием денежных средств имеет релейный характер

b( s) =

s < s„

bi, s > s0

(1)

Рассмотрим два близких момента времени / и / + Д/. В силу сделанных предположений изменение капитала фонда за время Д определяется следующим соотношением:

f-b(s) Д с вероятностью (1 - ХД ) + о( At),

As(t ) = <! (3)

[x - b( s)At с вероятностью ХДф( x)dx + о(Д ),

откуда, считая функцию p(s, t) дифференцируемой по s и t и используя формулу полной вероятности [3], получим [1], что p(s, t) удовлетворяет уравнению

л , ч ,, ч dp(s, t) dp(s, t)

Xp(s, t ) = b(s) +

ds dt

xj p(s - x, t)ф(x)dx. (4)

Решение уравнения (4) должно удовлетворять начальному условию

P(s,0) = Po (s),

(5)

где p0 (s) - заданная функция, и условию нормировки

j p(s, t )ds = 1.

(6)

для некоторого значения капитала 50. Так как фонд не имеет целью получение прибыли, то естественно считать, что при 0 <9 < 1

b0 = (1 - 9)Xa, b1 = (1 + 9)Xa.

(2)

Таким образом, при s < s0 фонд расходует в среднем меньше, чем собирает, а при s > s0 расходует в среднем больше денежных средств, чем в него поступает. Величина 9 имеет тот же смысл, что и нагрузка страховой премии в задачах страхования [2]. В дальнейшем считается, что 9<< 1.

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

Плотность распределения капитала фонда

Обозначим через p(s, t) плотность распределения вероятностей капитала фонда s в момент времени t .

Наконец, интегрируя уравнение (4) в пределах (-œ,+œ) и учитывая соотношения (1) и (6), получим, что при s = s0 должно выполняться условие

bop(s0 - 0, t) = bip(s0 + 0, t). (7)

В дальнейшем считаем, что s0 = 0, чего всегда можно добиться перенеся начало отсчета.

Пусть 9 << 1. Решение уравнения (4) будем иметь в виде

p(s, t ) = 9y(9s, X92t, 9).

(8)

Подставив выражение (8) в уравнение (4) после замены переменных

г = 95, т = Х92/, (9)

получим уравнение относительно функции у( 2, т, 9)

Ху(г.,,е> = bI zV1^zT'0)

0i dz

-X02 dy(z, t, 0) dr

- X J у (z - 0x, т, 0)ф( x) dx.

(10)

Пусть Е(2,ю) = |у(2,т)е ютйт - преобразование

0

Лапласа функции у(2, т). Применяя преобразование Лапласа к уравнениям (15) и (16), получим

Считая функцию у(2, т, 9) дважды дифференцируемой по 2 , можно записать

у(z - 0x, т, 0) = у(z, т, 0) - (z, т, 0)0x +

•• 02 x2

+ у (z, т, 0)-----------+ о(02).

z 2

(11)

a2 д

d

0 —2 У( z, т, 0) + а0 — у( z, т, 0) -

dz

2 ' dz2

-02 — у( z, т, 0) + o(02) = 0, дт

(12)

а2 d2F(z, ю) dF(z, ю)

2 dz2

+ а-

dz

--raF(z, ю) = 0, z > 0, (20)

a2 d2F(z, ю) dF(z, ю)

2 dz2

dz

-fflF(z, ю) = 0, z < 0. (21)

Подставляя разложение (11) в уравнение (10), получим, что в области 2 > 0

Рассмотрим область 2 > 0 . Характеристическое уравнение уравнения (20)

—к2 + ak -ю=0 2

имеет корни

кДю) =

-а + -.J а2 + 2а 2

к2(ю) =

-а --J а2 + 2а2'

-. (22)

а в области z < 0

а d2 d

—02 —- У( z, т, 0) - а02 — у( z, т, 0) -2 dz dz

-02 — у( z, т, 0) + o(02) = 0. dr

Пусть существует

у(z, т) = lim у (z, т, 0).

(13)

(14)

а2 d2у(z, т) dy(z, т) dy(z, т)

2 dz2

- + а

dz dr

а2 d2у(z, т) dy(z, т) _ dy(z, t)

2 dz2

dz dr

z > 0, (15)

z < 0. (16)

Поэтому

F (2, ю) = СДю)/1^ 2 + С2(ю)ек2( ю) 2.

Так как при 2 F(2, ю) ^ 0, то С1 (ю) = 0 . Та-

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

ким образом, при 2 > 0

F (z, ю) = C (ю)ек2( ю) z.

(23)

Переходя в уравнениях (12) и (13) к пределу при 1^0, получим

Рассмотрим область 2 < 0 . Характеристическое уравнение уравнения (21)

—ш2 - ат -ю = 0 2

имеет корни

а-+а2 + 2а 2 ю а + \/а2 + 2а,ю

ш1(ю) =-----*-—, ш2(ю) =----------------------------^-2. (24)

Условие нормировки (6) при этом перепишется в виде

| у(z, r)dz = 1,

0^о 0 0 I 0

(17)

условие сшивания (7) принимает вид

у(0 - 0, /) = у(0 + 0, /) (18)

и, наконец, начальное условие (5) дает

у(z,0) = lim1 Р0I| = S(z). (19)

Поэтому

F(2, ю) = й1(ю)е"1(и)2 + й2(ю)ет2(2 .

Так как при 2 ^ -оо F(2, ю) ^ 0, то й1 (ю) = 0 . Таким образом, при 2 < 0

F (z, ю) = d (ю)еш2(ю) z.

(25)

Далее условие сшивания (18) дает й(ю) = С(ю). Наконец , из условия нормировки (17) получим, учитывая, что к2 (ю) = -т2 (ю),

C (ю) =

m2 (ю) 2ю

Таким образом, окончательно

F ( z, ю) =

= m2 (ю) -т2(ю)|z|

(26)

Пусть у( z) = lim у( z, t).

получим, что

у( z) = — exp

а2

2—|z|

V —2 )

(27)

p(s,t ) =

_ (| s|+>—te)J 2Àa,t

yjJжkaJ,

2а2

2a0| s|

Erfc

фх—jt )

. (30)

Учитывая предельное соотношение между функцией у( 2, /) и ее преобразованием Лапласа F (2, ю) [4]:

lim у( z, t ) = lim oF ( z, ю),

t^œ z^0

Зная плотность распределения капитала фонда, можно определить такие его характеристики, как вероятность неплатежеспособности фонда и вероятность повышенных выплат. Рассмотрим только стационарный режим, когда плотность распределения капитала фонда определяется соотношением (28).

Фонд не может производить текущие выплаты, ко -гда капитал становится отрицательным (5 < -50 с учетом переноса начала отсчета). Из соотношения (29) вероятность этого события

Ро = p{S <_S0}= 2exP

(31)

Фонд производит повышенные выплаты, когда его капитал 5 > 0 (с учетом переноса начала отсчета). Из соотношения (29) вероятность этого события

При произвольных значениях / функция у(2, /) определяется как обратное преобразование Лапласа функции F(2, ю):

у(z,t) = —f F(z, o)eotdю . 2nj Ji

После несложных преобразований функцию F (г, ю) (26) можно переписать в виде

_—z

F(z,ю) =^L= e —2 -J— e-^

где

4)—) yfp _p

. a2 V2|z e

P =------, a = —f=^ , Р = ю + р

2a,

4—2

откуда [5]:

_ (| ^+—t )2

y(z,t) = , 1 e 2—2 +-^e “2 Erfc

■yj2na2t 2a2

2aN (і i Л

|z| _ at

У4J—J_t )

(28)

2 ^ 2

где Erfc(x) = —=• f e- dt. Vn x

. . 0a ( 2a s 0

p( s) =---------------------------------exp-L-1—

a2

V a2 )

+ o(0),

(29)

Pi

= p{s > о}=2.

(32)

Задавая вероятность p0, можно определить допустимые 9 и s0, определяющие заданную вероятность неплатежеспособности фонда.

Плотность распределения капитала фонда в стационарном режиме

В стационарном режиме полученное выше выражение (29) для плотности капитала фонда p(s) может быть уточнено. В соответствии с (4) в стационарном режиме плотность p(s) определяется уравнением

Xp(s) = b(s) dp(,y) +Х|p(s_ x)<p(x)dx . (33)

ds n

Решение уравнения (33) будем искать в виде

p(s) = 9y(9s, 9). (34)

Функция у ( z, 9), очевидно, удовлетворяет уравнению

Ху( z, 0) = 0b

z 1dy(z, 0) 0 ) dz

-Х|у(z _0x, 0)ф(x)dx. (35)

Переходя от функций у(z) и у(z, t) к плотностям распределения p(s) и p(s, t) соответственно, получим, учитывая (8) и (14),

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

dy(z, 0) + a2 d2у(z, 0) 0a3 d3y(z, 0)

dz 2 dz2 6 dz3

а в области z < 0

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

+ o(0) = 0,(36)

ду(г, 9) а2 д2у(г, 9) а3 д3у(г, 9)

дг

дг2

дг3

+ 0(9) = 0, (37)

где а3 = М {х3} . Пусть теперь

У(г, 9) = у 0(г) + у 1 (г)9 + 0(9).

(38)

Подставляя разложение (38) в уравнения (36) и (37) и приравнивая коэффициенты при степенях 9 , получим

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

ф(х) = -Хгехр [- х

Решение уравнения (33) в этом случае имеет вид Гё ехр^я), 5 < 0,

Р(?) = \ (46)

С1 ехр^я) + С2 exp( 125), 5 > 0,

+ а-

2 ёг2

а2 ё2у0(г) ёу0(г)

■ = 0, г > 0,

где

(39) К =

= 0, г < 0;

-(3 - 49) + У 9 - 89 4Ц1 -9)

-(3 + 49) + л/ 9 + 89

2 ёг2 ёг

а функция у1 (г) удовлетворяет уравнениям

а2 ё2у1(г) ёу1(г) а3 ёЗу0(г)

0! = ё а(1 -9)2(^ - ^-9 , С, = ё

4Х (1 + 9) а(1 -9)2 - 2-9

а(1 + 9)222 - г1)

- + а-

2 ёг2 ёг 6 ёг3

а2 ё2у1(г) ёу1(г) а3 ёЗу0(г)

-- а-

2 > 0, г < 0.

(40)

а(1 + 9)2( г1 - г2) и параметр

1 а(1 -9)2(К1 -г2)-9 а(1 -9)2(К1 -г1)-СП

ё =

К1 а(1 + 9)2(г1 -г2)г1 а(1 + 9)2(г2 -г1)г2

2 ёг2 ёг 6 ёг3

Условие нормировки (17) дает теперь

+да +да

| у0 (г)ёг = 1, | у1 (г)ёг = 0,

-да -да

а условие сшивания (7) перепишется в виде у с(0 - 0) = у0(0 + 0),

уД0 - 0)-у 0(0 - 0) = у1(0 + 0) + у„(0 + 0).

(41)

Графики функции р(5) (46) и ее аппроксимаций р1 (5) и р2 (5), полученных по формулам (29) и (45) соответственно, приведены на рис. 1. Параметр а = 1, параметр 9 = 0,1; 0,5.

(42)

Из соотношений (39), (41), (42) следует, что функция у0(г) определяется соотношением (27):

у0( г) = — ехр

ап

2а |г|

(43)

2

0=0.1

/

\ /

' Р1

? 4 / / N ч ■■■■■/

!'/ ■ V Р

/ / /

/ Р2

а решение системы (40), удовлетворяющее условиям (41), (42), имеет вид

уД г) =

V 3а24 4а3 а3 3а2

4а3 а а

3 -г--------

2

а Л

-I +-----

ехр

ехр

2

г > 0,

г < 0,

(44)

откуда при малых 9 плотность распределения капитала фонда

Р( 5) =

9а "

1 -9

а2 _ V

9а " (

1 + 9

а2 _ V

1 -

4а3а

2

-95

3а23

4а3а2 1 + —^ 95 3а2

V (

ехр

/_ V

Л" (

ехр

/_ V

2а9я

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

а

2а95

а2

5 > 0,

2

Л

(45)

5 < 0.

Рис. 1

Плотность распределения периода неплатежеспособности

Неплатежеспособность фонда наступает при 5 < 0. Обозначим через /(,?) продолжительность периода неплатежеспособности фонда при условии, что в начале периода капитал фонда равен 5 (естественно, 5 < 0). Пусть прошло времени Д/. За это время капитал фонда изменится на величину Д5 и

И 02 З2 f (v, z, 0)

+ Ха0

df (v, z, 0)

t (s) = At +1 (s + As).

(47)

(1 - 0)Ха дР(u, 5) + (и + X)Г(и, s) -

дs

- 5 да

-XI Г (и, 5 + х)ф( х)с1х -Х| ф( х)с1х = 0.

0 -5

Решение уравнения (48) будем искать в виде

Г (и, 5) = / {и-, 05, 01 .

(48)

(49)

Подставляя выражение (49) в уравнение (48) после

замены переменных г = 9-5 , V = и нение относительно / (V, г, 9)

2 dz2 dz

-v02 f (v, z, 0) + o(02) = 0.

Пусть существует

f (v, z) = lim f (v, z, 0).

(52)

Переходя в (51) к пределу при 9^0, получим, что функция / (V, г) определяется уравнением

Обозначим через Е(и, 5) = М {е-и>(5)} характеристическую функцию величины /(5). Из соотношения (47) имеем, очевидно, Е (и, 5) = е-“Д'МД5 {Е(и, 5 + Д-)} .

Или, учитывая (3), при малом Д/ и 5 < 0

F(u, s) = e uAt x

-s да

(1 -XAf)F(u, s - b0At) +XA I F(u, s + х)ф(х) dx+XA | ф(х)ск

0 -s

+ o(At),

так как при 5 > 0 /(5) = 0 . Считая функцию Е(и, 5) дифференцируемой на 5 и переходя к пределу при Д/ ^ 0 , получим уравнение, определяющее Е(и, 5),

И 52 f (v, z) +Xa df (v, z)

2 5z2

dz

-- vf (v, z) = 0 (53)

с вытекающим из (50) начальным условием / (т,0) = 1. Характеристическое уравнение уравнения (53)

Ха2

к + Хак - v = 0

имеет корни

-Ха + л/х2а2 + 2Ха7 V -Ха - л/х2а2 + 2Ха,у

К =--------*------------—, К =----------^------------2 .(54)

Ха2 Ха2

Учитывая, что |Е(и, 5) < 1 и, следовательно, \/(V, 5)| < 1, получаем, что решение уравнения (53) имеет вид

f (v, z) = exp [(v) z ],

(55)

откуда при 9 << 1 характеристическая функция Е(и, 5) определяется соотношением

F (u, s) = exp

кі

0s

получим урав-

Плотность распределения капитала 5 при условии, что 5 < 0, как вытекает из соотношений (29) (31), в стационарном режиме имеет вид

9(1 - 9)Ха д(^ г 9) + (9^ + Х)/(V, г, 9) -дг

-У9 да (50)

-Х | /(V, г + 9х, 9)ф(х)ёх - Х | ф(х)ёх = 0.

0 -!9

Будем предполагать, что при г >> 1 ф(г) < ^

для некоторых ц и у, а функция /(V, г, 9) дважды дифференцируема по г. Тогда, раскладывая / (V, г + 9х, 9) в ряд Тейлора и ограничившись первыми тремя членами разложения, получим

p(|s < 0) =-----

Поэтому безусловная характеристическая функция периода неплатежеспособности

F(u) = | F(u,s)p(|s < 0)s = —

-да 1 -

+ -\/Г + au

(56)

где а = 2а2 Ха292 .

Вычисляя обратное преобразование Лапласа от функции (56), получим [5], что при 9<<1 плотность распределения периода неплатежеспособности имеет вид

2

e

2

р(/) = .2 ехр |—-|-—Бйе л/па/ V а У а

а

V У

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

1

ф(х) =- е-ш .

а

Плотность распределения периода неплатежеспособности определяется при этом соотношением [1]:

Х 9 да 1

р0 (/) = ^= Г- /1 (2^ТГ-0х)е^Х(2-9)хёх, (58)

•у 1 -9 ' х

где /1( х) - модифицированная функция Бесселя первого порядка. Графики функции р0 (/) (58) и ее аппроксимации р(/) (57) при Х = 1 и 9 = 0,1; 0,5 приведены на рис. 2.

у Р0,©=С .5

\ \ Р0, Р1 ^0=1 Р1 ,©=0 5

-

/(5) = Д/ + /(5 +Д5) .

(59)

или, учитывая (3), при малом Д/ в области 5 > 50 получим

Е(и, 5) = е~“Д' X

да

(1 - ХД/) Е (и, 5 - ЬД) + ХД/Г Е (и, 5 - Ь1Д/ + х) ф( х) ёх

0

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

+о(Д/).

Переходя к пределу при Д/ ^ 0 и считая функцию Е(и, 5) дифференцируемой по 5 , получим уравнение, определяющее Е(и, 5):

(1 + 9) Ха дЕ (и 5) + (и + Х) Е (и, 5) -

д5

да

-Х Г Е (и, 5 + х)ф( х)ёх = 0.

(60)

При 5 = 50 величина 50 - Ь1Д/ < 0 . Поэтому при

Е(и, 5) = е~“Д' X

да

(1 - ХД/) + ХД/Г Е (и, 5 + х) ф( х) ёх

+ о( Д/).

Переходя к пределу при Д/ ^ 0 получим отсюда граничное условие в точке 5 = 50

Е(и, 50) = 1. (61)

Решение уравнения (61) будем искать в виде

Е(и, 5) = /\ —, 95

Рис. 2

Плотность распределения периода повышенных выплат

Период повышенных выплат наступает, когда капитал фонда 5 > 50. Обозначим через /(5) продолжительность периода повышенных выплат, если в начале периода капитал фонда равен 5 (5 > 50). Пусть прошло время Д/. За это время капитал фонда изменится на величину Д5 и

(62)

Подставляя выражение (62) в уравнение (60) после замены переменных г =9-5, V = и |92, получим уравнение относительно / (V, г, 9)

(1 + 9)9Ха д/(^ г 9) + (9^ + Х)/(V, г, 9) -дг

да

-Х Г / (V, г + 9х, 9)ф( х)ёх = 0.

(63)

Обозначим через Е(и, 5) = М {е и(5)} характеристическую функцию величины /(5). Из соотношения (59) имеем

Е(и, 5) = е^М^ {Е(и, 5 + Д5)} ,

Считая функцию /(V, г, 9) дважды дифференцируемой по г , раскладывая /(V, г + 9х, 9) в ряд Тейлора и ограничиваясь первыми тремя членами разложения, получим

Ха2 92 д2/(V, г, 9)-Ха92 д/(V, г, 9) 2 дг2 дг

^92 / (V, г, 9) + о(92) = 0.

(64)

Пусть существует

f (v, z) = lim f (v, z, 0).

(65)

Переходя в (64) к пределу при 9^0, получим, что функция / (V, г) определяется уравнением

Ха^ _хй -V/ (V, г) = 0. (66)

2 дг дг

Характеристическое уравнение уравнения (66)

Ха2

к - Хак - v = 0

имеет корни

Ха -л/Х2 а2 + 2Ха^ Ха + л/Х2 а2 + 2Ха7v

к = *__________^ К =__________________—.

Ха2 Ха2

Учитывая, что |/(V , г)| < 1, получим, что решение уравнения (66) имеет вид

/(V, г) = СМехр (к^) г) ,

откуда при 9<<1 характеристическая функция Е(и, 5) определяется соотношением

F (и, 5) = Cl-Url exp [к! i-Ur 105 I.

Или, с учетом условия (61),

F (и, 5) = expl кЛ — 10(5 - 5о) I .

Плотность распределения капитала 5 при условии, что 5 > 50, как вытекает из соотношений (29) и (33), имеет вид

I ч 2а0

Р(5 5 > 5о) =-------------e

а2

2а0(5 - 5о)

Поэтому безусловная характеристическая функция периода повышенных выплат

F(и) = | F(и, 5)Р (5 > 50 )5 = — 51

+ >/Г + аи

(68)

где а = 2а2 Ха292 .

Вычисляя обратное преобразование Лапласа от функции (68), получим, что при 9<< 1 плотность распределения периода повышенных выплат имеет вид

2 [ t | 2

p(t) = -------exp l------I------Erfc

л/nat v aJ а

( . 1Л

112

(69)

Таким образом, в рассмотренном выше симметричном случае, определяемом соотношением (2), при 9 << 1 плотности распределения периодов повышенных выплат и неплатежеспособности совпадают.

а

ЛИТЕРАТУРА

1. Лившиц К.И., Шифердекер И.Ю. Математическая модель деятельности некоммерческого фонда при релейном управлении капиталом //

Вестник Томского государственного университета. Приложение. 2006. № 18. С. 302-308.

2. Глухова Е.В., Змеев О.А., Лившиц К.И. Математические модели страхования. Томск: Изд-во Том. ун-та, 2004. 180 с.

3. ГнеденкоБ.В. Курс теории вероятностей. М.: Наука, 1969. 400 с.

4. ЛаврентьевМ.А., Шабат Б.В. Методы теории функций комплексного переменного. М.: Наука, 1973. 736 с.

5. Бейтмен Г., Эрдейи А. Таблицы интегральных преобразований. М.: Наука, 1969. Т. I. 344 с.

Статья представлена кафедрой прикладной математики факультета прикладной математики и кибернетики Томского государственного университета, поступила в научную редакцию «Кибернетика» 30 мая 2006 г.

i Надоели баннеры? Вы всегда можете отключить рекламу.