CC BY
31
^ Научно-технические ведомости СПбГПУ. Физико-математические науки № 1(189) 2014 ^ УДК 517.968.22
А.В. Кондрашков ', Ю.А. Пичугин 2
' Санкт-Петербургский государственный экономический университет 2 Российский государственный педагогический университет им. А.И. Герцена
идентификация и статистическая проверка устойчивости
модели вольтерры
В статье рассматриваются некоторые аспекты модели Вольтерры. Основное внимание уделяется определению области устойчивости в пространстве параметров модели и статистической проверке устойчивости.
модель вольтерры, идентификация, устойчивость, статистическая проверка.
Введение
В задаче прогнозирования того или иного процесса (экономического, социального, биологического и т. п.) единственным источником фактической информации служат данные наблюдений, представимые в виде таблиц чисел. Здесь прослеживается связь с задачами анализа экспериментальных данных. Чтобы преодолеть ситуацию «черного ящика», в ряде случаев бывает целесообразно выбрать какую-нибудь известную математическую модель в качестве отправного пункта поиска подходящей модели, требуемой для интерпретации наблюдаемых данных. Например, в теории экономических циклов имеется динамическая модель Гудвина [1] на основе системы Лотки — Вольтерры [2], связывающая пару величин. Однако модель такого вида является слишком «жесткой», поскольку демонстрирует единственный сценарий (замкнутые траектории). Простое продолжение на основе общей системы Вольтерры ведет к более вариативной модели [3]. При этом важная величина, характеризующая эластичность в паре переменных, остается в пределах класса дробно-линейных выражений. Здесь следует подчеркнуть, что выбор типа модели всегда является прерогативой предметного исследователя и не требует особых усилий в плане обоснования, так как ошибка выбора непременно приводит к ситуации, когда дисперсии оценок параметров относительно велики, а сами оценки параметров не проходят статистического тестирования, что означает отклонение модели.
В процессе реализации выбранного типа модели возникает весьма обширный круг задач. Однако ниже мы ограничимся рассмотрением лишь следующих этапов: оценка параметров модели; определение (в пространстве параметров) области устойчивости точного решения; статистическая проверка принадлежности параметров области устойчивости. Остальные вопросы предполагается рассмотреть в дальнейших публикациях.
Постановка задачи
Предположим, что для пары неких предметных величин х и у требуется выяснить функциональную зависимость между ними на основании исходных данных, представи-мых в виде пары временных рядов наблюдений £ и п. Пусть известно, что данные наблюдений имеют шаг 5 по шкале времени t на промежутке t0 = 0 < t < ^. Временные ряды Е = (Е,) и п = ("Л,-) можно рассматривать как сеточные функции на равномерной сетке
ю5 : = {t¡ = ,5 : , = 0, т}.
При дальнейшей обработке данных нам понадобятся следующие числовые характеристики:
__т __т
Е := (1/ Л := (1 /
¡=1 ¡=1
т ___
соу(Е, Л) : = (1 / т)£ (Е,- - Е)(Л,- - л)-
,=1
С целью обнаружения и обоснования функциональной зависимости широко
применяют методы регрессионного анализа. Традиционные подходы связаны с поиском уравнений регрессии на основе явных форм:
у = У (х; О или х = X (у; О,
где С — списки определяемых параметров.
Между тем в каких-то случаях путь рассеяния данных может свидетельствовать в пользу выбора неявной формы:
Ф(х, у; о = о.
Обобщенный подход связан с описанием на основе параметризации. С этой точки зрения графиком искомой функциональной зависимости будет некая геометрическая кривая, представимая как множество всех эквивалентных путей
{Р ^ (х(р), У(р))},
параметризующих эту кривую.
Достаточно гладкая общая кривая, заданная уравнением к(х, у) = 0, допускает параметризацию с помощью решения задачи Коши для «почти гамильтоновой» системы вида
(х , ,дк, ч — = -а(х, у) —(х, у); аТ ду
= -а(х, у) д^ (х, у). (Т дх
Начальные условия при Т = 0 имеют вид:
х(0) = хо, у(0) = уо; Л(хо, уо) = о.
С другой стороны, параметризованная кривая, заданная уравнениями (х, у) = (х(Т), у(Т)), имеет представление в виде локально конечного объединения общих кривых, уравнения которых могут быть получены путем исключения параметра Т из параметрических уравнений.
Отдельно взятый путь
Т ^ (х(Т), у(Т)),
где Т — параметр времени, может служить описанием процесса [4, 5].
Динамическая модель процесса
Пусть через Т обозначено время. Пусть
(1)
гладкий путь (х(Т), у(Т)) служит изображением наблюдаемого процесса. Тогда можно принять динамическую модель (с непрерывной шкалой времени) на основе системы обыкновенных дифференциальных уравнений (ОДУ):
\х = ^(х, у);
I у = Щ х, у),
где F = (^1, — гладкое векторное поле.
Как было отмечено выше, выбор выражений ¥1 и ^2 остается за предметным исследователем в соответствии с предпочтениями предметной области. Например, нередко авторы работ по математическому моделированию экономических процессов [4] идут по пути заимствования подходящих динамических моделей из биологии. Фактически эти заимствования с учетом тех или иных традиционных предположений [1] определяют общий вид ожидаемых динамических систем.
Положения равновесия (они же — стационарные точки) определяются как решения векторного уравнения F (х, у) = о.
Изоклина векторного поля F — это полезная линия, вдоль которой выполняется условие вида
V) = 0,
т. е. Б ± V при Б ф 0, где V ф о — фиксированный вектор.
Ниже используются Л-изоклины («горизонтальные», когда ^2 = о) и у-изоклины («вертикальные», когда F1 = о).
Фактически любая стационарная точка системы ОДУ (1) реализуется как точка пересечения у-изоклины и Л-изоклины.
Специфика динамической модели. При моделировании процессов в ряде предметных областей (к числу которых относятся экономика, социология, биология, химия, и т. п.) встречаются неотрицательные предметные величины, которые могут принимать только неотрицательные значения. Предполагается, что х > о и у > о — именно такие величины. Тогда можно несколько уточнить общий вид уравнений динамической системы (1), чтобы автоматически иметь описание в положительном конусе:
:= {(х, у): х > о, у > о}.
С этой целью достаточно принять
X = х/х(X, у), у = у/2(X, у). (2)
Заметим, что сделанный выбор соответствует специфике популяционной динамики по Колмогорову [3].
На основании данной динамической модели можно написать ОДУ:
х^у = /¡(х, у) у йх /2(х,у)'
Выражение левой части этого уравнения обычно называют эластичностью. Таким образом, эластичность выражается формулой:
Э (X, у) =
/2( х, у) /¡( X, у)
Идентификация модели. Ниже используется описание в положительном конусе И+. В таком случае динамическая модель (2) приобретает следующий вид:
- = /¡(X, у);
X
- = /»(X, у).
I у
Очередной этап сводится к нахождению правых частей этих дифференциальных уравнений по имеющимся наблюдаемым данным, заданным в виде пары временных рядов £ и п.
При отсутствии каких-либо иных соображений можно предположить, что предметный исследователь сделал выбор, приняв решение о поиске /¡( X, у) и /2( X, у) в виде многочленов от (X, у) степени < к.
В линейном случае (к = 1) порядок действий состоит в следующем. Наряду с имеющимися наблюдаемыми временными рядами £ и п можно ввести для величин X / X и у / у вычисляемые временные ряды Ф = (ф;) и V = (¥;):
Ф :=
П.
П
Временные ряды £ и п рассматриваются как сеточные функции, так что £. и п; — это стандартные обозначения разностных производных:
(?;), ; (п) = .
Далее, с целью определения коэффициентов многочленов правой части можно применить метод наименьших квадратов (МНК) по схеме «2 * 3». В результате получим дифференциальные равенства вида:
X у
- = а + риX + Р12у; - = а» + Р21X + Р22у. (3)
X у
Для определения двух троек чисел (а1, Р11, Р12) и (а2, Р21, Р22) стандартные оценки по методу МНК в матричной форме имеют вид:
В = БС-1; (а1, а2)т = (ф, - В(ф, (4)
С использованием сокращенных обозначений (с?п := соу(?, п)) матрицы в (4) вводятся следующим образом:
(в
В=
11
УР21 в2
Гс
12 ; С := ??
22 У с V П?
' с. с Л
= Ф? ФП
с ? с V ¥П У
пп у
В итоге соотношения (3) могут быть записаны в виде системы ОДУ:
(5)
[X = X (а1 + Ри X + Р12 у); [у = у(а2 + ^ + Р22 у).
Важно отметить, что в популяционной динамике система ОДУ такого вида известна как общая система Вольтерры (ОСВ).
Выбранная динамическая модель (5) характеризуется дробно-линейной эластичностью:
э (X, у) =
а2 + Р21X + Р22 у
а1 + Pих + Р12 у
Расширенная матрица параметров может быть записана следующим образом:
В := ГР10 := а1 Рп Р1
Л
ув2о : а2 Р21 Р22 У
Вычеркивая столбцы этой матрицы, вычисляем миноры второго порядка:
йо := ае1 В = Р11Р22 -Р12Р21; й1 := а1Р 22 — а2P12, й2 := а1Р21 — а2Р11.
Полученная динамическая модель изображается «точкой» в пространстве матриц В размера 2 * 3. Равенства (4) позволяют перейти к укороченному описанию, так что модель процесса однозначно определяется квадратной матрицей В размера 2 * 2.
Чтобы развернуть процесс во времени, следует найти решение (х(Т), у(Т)) задачи Коши для системы ОДУ (5) при Т > о с использованием начальных условий:
х(о) = £о, у(о) = По.
(6)
Влияние матрицы параметров модели на характер процесса
Положения равновесия и изоклины. По
итогам проделанной идентификации векторное поле F приобретает компоненты:
Fl := х(а + рих + р12у);
¥2 := у(а,2 + Р21 х + Р22у).
Для более наглядного представления положений равновесия можно выделить полезные линии:
у-изоклины Оу и Ь^: а1 + рих + Р12у = о; Л-изоклины Ох и Ь2: а2 + Р21 х + Р22у = о. Тривиальные стационарные точки расположены на координатных осях Ох и Оу:
Ро := (о, о); Р1 := (х1,о) е Ь п Ох;
Р2 : = (о, у2) е Ь, п Оу.
По умолчанию подразумевается, что Ь ф Ох и Ь2 ф Оу. Условия реализации точек Р1 или Р2 :
если Р11 ф 0, то Ь п Ох = {Р1} и
х1 = —а1 / Рц;
если Р22 ф 0, то Ь, п Оу = {Р2} и
у2 = -а2 / Р22.
Нетривиальная стационарная точка
Р* := (х*, у») е Ь1 п Ь
определяется в случае Ь ф Ь> как решение системы линейных уравнений:
|а + ри х* + Р12 у, = о; [а2 + Р21 х* + Р22 у* = о. Условие реализации точки Р* (существование и единственность) состоит в невырожденности матрицы В:
Ь п Ь = {Р*} о ёе! В ф о.
(7)
По формулам Крамера имеем:
(о х* = -(1, (о у* = (2.
Процесс и его модель называются невырожденными, если ёе! В ф о. В настоящей работе не обсуждаются иные возможности, ведущие к возникновению стационарных решений.
Центрированные процессы. Процесс и его модель называется центрированными, если Ь1п Ь2 = {Р*} и Р* е И+ (точка Р* является внутренней точкой положительного конуса.)
Условие реализации центрированного процесса:
(1(о < 0, (2(о > о.
Остальные процессы. Нецентрирован-ные случаи появляются всякий раз, когда имеет место одно из двух условий:
либо (1) ае1 В ф 0, но Р* г И+; либо (2) ёе! В = о.
Генезис модели Вольтерры. Невырожденная модель М с матрицей В ведет свое происхождение от некой невырожденной модели М о с матрицей Во, если их матрицы принадлежат связной открытой области 0 пространства матриц 2 * 2, а соответствующие стационарные точки ( Р* и Ро,* ) принадлежат заранее заданному выпуклому открытому конусу К на плоскости переменных (х, у). При этом в области 0 не должно быть вырожденных матриц. При рассмотрении центрированных моделей следует принять К := И+.
Пусть Мо — это какая-либо известная модель (например, имеющая уникальное наименование и зафиксированная в литературных источниках). Если область 0 — максимальная среди всех, отвечающих указанным требованиям, то такая область 0 интерпретируется как клан модели Мо.
Расслоение пространства параметров модели. Этот подход способствует визуализации и позволяет выявить детально структуру границы того или иного клана.
Сначала введем структурированное арифметическое пространство ИВ как пространство всех квадратных матриц В размера 2 * 2. Возьмем за основу «перекрестное» представление:
И 4 ~ И 2 у И 2
В 11,22 12,21
где И2и2 и И22 21 состоят из пар р = (Ри, Р22)
и 9 ^Р^ Р21).'
Тем самым мы можем рассматривать И В как объединение всех слоев
^ := КП,22 * {9}.
В силу данного расслоения пространства ИВ возникает расслоение какой-либо области и в виде семейства слоев
ид := и п!?(д е рг и), где рги — это проекция и в пространстве
12,21
При этом возникает сопутствующая задача: требуется определить край ди? каждого слоя ич при всех д е рг и.
Область и и ее граница ди могут быть восстановлены путем объединения всех слоев области и.
Две области центрированных моделей. В случае невырожденной модели из сравнения формул (4) и (7) находим:
(х*, у*)г = (I,П)г - В-Чф, ;
х =| в22 Ф в12 V . у =~ в11 У-Р21Ф х* ёе! В ' у* = П ёе! В '
Чтобы несколько упростить дальнейшие выкладки, введем многочлены второго порядка от р..:
и* =| ёе!В - (Р22Ф-Р12 V);
у* = п ёе! В - (рпу-р21 Ф).
Модель с матрицей В реализуется как центрированная, если (и только если) выполняется один из двух вариантов условий:
1) когда ёе! В > о и тт(и*, у*) > о;
2) когда ёе!В < о и тах(и*, у*) < о.
Пусть область и' пространства ИВ определяется условиями варианта 1. Ее граница ди' может быть составлена из кусков трех алгебраических гиперповерхностей второго порядка:
с1е1 В = 0; и* = (3, у* = о.
Если же придерживаться расслоения, то в проекции на пространство И2122 край
дид произвольно взятого слоя и' формируется из кусков ветвей трех гипербол:
с1е1 В = 0; и* = 0; у* = о.
Особую роль играют места соединения разных кусков границы. В частности, чтобы определить точки пересечения гиперповерхностей и* = о и у* = (3, можно составить систему:
Г| ёе! В + Р12 У = р22 ф; [п ёе! В + Р21Ф = Ри V. Как следствие этой системы выводится уравнение сборки границы:
ёе!(В)(|п ёе! В + Р12 vn + Р21 ф| - ФV) = о. В дальнейшем выражение
:= Р12 Vп + P2l Ф|-ФV \п
рассматривается как характеристика сборки границы (из кусков).
Все эти рассуждения переносятся без изменений на случай другой области и", определяемой по варианту 2.
Рамочные условия. Дальнейшие выкладки будут правомерны, если исходные временные ряды £ и п отвечают определенным требованиям.
Чтобы обеспечить специфику, по умолчанию подразумевается, что указанные временные ряды — положительные: 11 > о и п. > о для всех 1 = 0, т. Таким образом, всегда имеем |> о и п> о. _ _
Однако средние значения ф и V вычисляемых рядов ф и у могут повлиять на схему сборки границы области центрированных моделей (в каждом из двух случаев). Выделим условия, руководствуясь следующим определением:
исходные данные (временные ряды £ и п) находятся в общем положении, если ФФ о и V ф о. _ _
Другие случаи, когда ф = о или V = 0> не рассматриваются в настоящей работе. В ситуации общего положения комбинация знаков ^п ф, V) может указывать на выбор первого из четырех главных кланов. Для каждого из этих кланов его образ при проекции в И2221 совпадает с первой из четырех компонент связности области
Р12Р21 ф 0. Ниже приводятся расчеты для одного из реальных кланов.
Клан модели Лотки — Вольтерры. В динамической модели (5) система Лотки — Вольтерры (СЛВ) реализуется, если
Р12Р21 < 0 и рп = Р22 = 0.
Имеются две возможности: либо q е (¿(а), либо q е (2т. Обозначения:
((а) : = к = (Р12, Р21): Р12 < 0, Р21 > 0},
6(4) := —((а).
Варианты (а) и (4) различаются лишь «направлением обхода» траекторий. Поэтому имеет смысл ограничиться рассмотрением одного из вариантов СЛВ, например (а).
Пусть для определенности в исходной ситуации общего положения имеют место следующие условия: ф > (3, V < 0. Тогда для всех q е ((а) при Р11 = Р22 = 0 выполняются неравенства:
и* = ?к + р 12 V > 0, V* = г|к + р21 ф > 0,
где к := —Р12Р21 > 0.
Тем самым все реализации СЛВ (а) — центрированные.
Все модели клана СЛВ (а), который далее обозначается через О, относятся к области ёе! В > 0. Характеристика сборки границы: х < 0.. Основная часть границы дО клана О формируется из трех кусков: Н1 := {ёе! В = 0, р11 > 0}; Н2 : = {V* = 0^ В > 0}; Н3 : = {и* = 0^ В > 0}.
Куски Н2 и Н3 соединяются вдоль многообразия Ь при Р11 < 0. Многообразие Ь размерности 2 определяется системой линейных уравнений:
|Р22 Ф- Р12 V =
[РпV-Р21Ф = 0.
Очевидно, что ёе! В = 0 на Ь. Если рассматривать произвольно взятый слой Од, то его край
дОq = (Hl)q ^ (H2)q ^ ^ .
В проекции на К2122 возникают следующие отождествления: (Н^ — ветвь гиперболы ёе! В = 0 (Р11 > 0); (H2)q — кусок
ветви гиперболы V* = 0(ёе!В > 0); (() — кусок ветви гиперболы и* = 0(с1е1 В > 0).
Пересечение Ь nEq проектируется в точку
Рь ) := (Р21Ф / V, Р12 V / Ф).
Остаточная часть границы состоит из кусков, допускающих расслоение при Р12Р21 = точнее, при q е д({ау
Условия устойчивости модели Вольтерры
Линеаризация в окрестности нетривиального положения равновесия. Перейдем к новым переменным:
X = X* +Ср у = у* +С2.
Система (5) преобразуется в систему обыкновенных дифференциальных уравнений:
К1 = ХфцС! +Р12^2) + с1(рЦс1 +Р12^2),
[С 2 = у*(Р21^1 + Р22^2) + С2(Р21^1 + Р22^2).
(8)
Для использования записей в матричной форме вводятся обозначения:
г := (Уи^ := X*PU, := у*Р2
С = (С1, ^)т.
Матричная запись системы (8) имеет вид:
С = ГС + С1ЗЕ(С1, С2)В с. (9)
Отбрасывая слагаемые вида 0(| £ |2), приходим к линеаризованной системе: £ = Г2. Инварианты матрицы Г выражаются следующим образом:
!г Г = Уп + У22 = ^Рп + У*Р22; ёе!Г = упу22 — 712У21 = X*у* ёе! Г.
Собственные значения матрицы Г могут быть получены как корни ее характеристического многочлена
ёе^Г — АЕ) = А,2 — А 1:г Г + ёе! Г.
Линеаризация в окрестности тривиального положения равновесия. Сначала вычислим для векторного поля Е его матрицу Якоби в произвольно взятой точке:
(р)т дх д^/ ду Л
^д/2 / дх д^2 / ду
(
Л
/1 + Р11х Р12х ^
ур21 у /2 + в22 у, Дальнейшая линеаризация системы (5) в окрестности каждой из точек Рк (к = 1, 2, 3) по известной схеме ведет к получению линейной системы:
= Гк г; К : = (ЩРк))
При возникновении необходимости все элементы матриц Гк могут быть выражены через элементы расширенной матрицы В.
Устойчивость и неустойчивость. Руководствуясь теорией Ляпунова [6], можно выделить достаточные условия устойчивости и неустойчивости положения равновесия на основании спектральных свойств матрицы линеаризованной системы. Применимость теорем Ляпунова в случае системы (5) обеспечивается тем, что отбрасываемые при линеаризации слагаемые являются однородными многочленами второй степени.
Какая-либо квадратная матрица А называется устойчивой, если ее спектр содержится в левой полуплоскости: Я.е X < о. Для матрицы размера 2 * 2 условие устойчивости равносильно выполнению двух неравенств: с^ А > 0, 1х А > о.
Нетривиальное положение равновесия характеризуется выполнением одного из трех случаев:
если ёе! Г> о и Г< о (матрица Г — устойчивая), то стационарная точка Р* асимптотически устойчива (локально);
если ёе!Г> о и !гГ> о (матрица (—Г) — устойчивая), то стационарная точка Р* вполне неустойчива (локально);
если с^ Г < 0, то стационарная точка Р* неустойчива (коль скоро она является седловой).
Во всех остальных случаях (когда ёе! Г = о или Г = о) требуются индивидуальные рассмотрения. В частности, если ёе! Г > о и с1е1 Г = 0, то выпуклая окрестность точки Р* заполнена замкнутыми траекториями [7], так что имеет место простая устойчивость.
Как для точки Р* , для каждой из точек Рк (к = 1, 2, 3) можно составить аналогичные условия устойчивости и неустойчивости. Формулировки условий здесь не пред-
ставлены, поскольку в рассматриваемых ниже задачах точки Рк себя не проявляют.
Области устойчивости и неустойчивости модели Вольтерры. Здесь и далее толкование термина «устойчивость» связывается с выполнением условий устойчивости нетривиального положения равновесия. В порядке очередной задачи требуется определить в пространстве ИВ связную открытую область О, состоящую из таких матриц В, для которых матрица Г — устойчивая.
Матрицы В и Г связаны соотношением Г = х*, у*)В.
Все дальнейшие расчеты выполнены в клане О, т. е. в клане СЛВ (а). Ввиду того, что
ёе! Г = х*у* ёе! В,
в клане О имеем: ёе! Г > о. В качестве области устойчивости можно принять область
О := {В е О : !г Г < о}.
Чтобы существенно упростить выкладки, введем выражение (многочлен третьей степени от р.):
р := ёе^Г) !г(Г) = и*Р11 + у*Р22.
С использованием р возникает другое представление области устойчивости: О := {В е О : р < о}.
Основная часть границы дО. области устойчивости О формируется из трех кусков по списку {Но, Н-, Н2}, где
Но : = {р = о^Г > о};
Н1- := Н1 п {р < о}.
После несложных преобразований можно написать:
р = ёе!(В)[рп| + р22П- (Ф + У)] +
+ к(Ф + V) + (РИР12 V + Р22Р21 Ф).
Пересечение гиперповерхностей
ёе!В = о и р = о имеет вид Хи N где многообразие N определяется условиями Р22 = Р12, Р11 = Р21. Куски границы Но и Н2 соединяются вдоль многообразия Ь при Р11 < о. Куски Но и Н1- соединяются вдоль многообразия N при Р11 > о.
Для произвольно взятого слоя Од области устойчивости его край имеет представление
Математика -►
Области устойчивости (О) и неустойчивости (50 и 51 и 52) для клана модели Лотки — Вольтерры в случае выполнения неравенств: р12 < О, Р21 > О
= (Яо) г и(Н-) ? и(Н2)д.
В проекции на 22 возникает еще одно отождествление, согласно которому (Я0)г — кусок алгебраической кривой (третьего порядка): р = 0 при ёе!В > 0. Пересечение N пЕ? проектируется в точку
Рм (Я) := (в21, в12).
Остаточная часть границы области устойчивости состоит из двух кусков, расположенных на гиперплоскостях Р12 = 0 и
Р2
0.
На рисунке показаны геометрические объекты в сечении слоем (индекс q опущен): О — область асимптотической устойчивости клана (?; 50 — область вполне неустойчивости клана (?; 51 и 52 — области, не задевающие (?, составленные из множеств, в каждом из которых имеет место одно из двух: либо модель нецентриро-ванная, либо она центрированная, но точка Р* — неустойчивая (когда положение равновесия является седловой точкой).
Устойчивость с точки зрения проверки статистических гипотез
В случае, когда полученные значения коэффициентов указывают на устойчивость процесса, описываемого системой (3), исследователю хотелось бы выяснить, насколько надежен результат, что, естественно, переводит задачу в область проверки статистических гипотез. До этого момента даже характерные для статистики формулы использовались лишь в смысле МНК, который, как известно, весьма популярен в области прикладной статистики.
Теперь нам следует подчеркнуть, что вычисленные значения параметров (коэффициентов) системы (3) являются оценками, которые, в принципе, носят случайный характер, и поэтому мы будем их обозначать как Р(>., а обозначения без «крышек» использовать для истинных значений, которые, естественно, считаем неизвестными. Будем предполагать, что оценки коэффициентов подчиняются нормальному
распределению. Поэтому ближайшая наша задача — это выяснение параметров этого распределения, а именно — матрицы дисперсий и совместных ковариаций.
Реализации левых частей системы (3), как и ранее, обозначаем через ф. и у., соответственно, а реализации переменных х и у — соответственно через £. и п., где , = 0, т. Перейдем к центрированным значениям Ф, = Ф, - Ф, где ф - среднее арифметическое. С остальными переменными поступим аналогично. Построим выборочную
матрицу левых частей системы (3):
(_ _ ^
А =
Фт
М7 т
Ф1 Ф2
Пусть Р — ортогональная матрица, такая, что РТААТ Р = Л, где Л - диагональная матрица. Тогда в строках матрицы Б = РТ А стоят реализации независимых, нормально распределенных величин, которые мы обозначим через и £2. Рассмотрим систему регрессионных уравнений:
К = ¿цх + ¿12 у + 61;
= ¿21 х + Ь22 у + е2.
Пусть через В7 обозначается (2 * 2)-матрица параметров системы (1о). Тогда МНК-оценка этих параметров следует выражению [8, 9]:
В7 = БХТ (ХХТ )-1,
где X - выборочная матрица центрированных реализаций переменных х и у, т. е.
(;
X =
Л
11 12 • 1т ЧП1 112 ••• Пт У
Остатки регрессии (1о) е.. (. = 1, 2) имеют нулевое среднее значение, а несмещенные оценки дисперсий этих остатков следуют выражению:
а2 = (т - 3)-1 £ (I, - ¿..11 - )2. (11)
I=1
Обращая внимание на то, что матрица коэффициентов при переменных в правой части (3) может быть получена по формуле
В = АХТ (ХХТ )-1, приходим к соотношению
В = РВ.
Р11 о Р12 о
о Р11 о Р12
Р21 о Р22 о
о Р21 о Р22,
Отсюда видно, что векторы
РТ = (Р11, Р12, Р21, Р22), ЬТ = (¿11, ¿12, ¿21, ¿22) связаны следующим соотношением [8, 9]:
Р = 0Ь, (12)
где ортогональная матрица () составлена из элементов Р, т. е.
е =
Здесь и далее для какой-либо матрицы В = (Р.) размера 2 * 2 через р обозначается ее развертка в столбец, так что
РТ := (ри, Р12, Р21, Р22).
Матрица дисперсий и взаимных кова-риаций вектора ¿ в силу независимости z1 и £2 имеет блочную структуру и выражается как
(д2(ХХТ )-1 О ^ ¿¿ О д2(ХХТ)-1 у
где каждый блок относится к паре коэффициентов соответствующего уравнения.
Из соотношения (12) следует, что матрица дисперсий и взаимных ковариаций вектора р равна
Т
V- =
В рассмотренной выше ситуации мы имеем устойчивость, когда точка р с координатами (Р11, Р12, Р21, Р22) лежит в области О, т. е. когда р е О. Вектор р является несмещенной оценкой, следовательно, этот вектор подчиняется нормальному распределению, т. е. р ~ N(р, У-). Функция плот-
ности этого формулой
распределения выражается
/(Р) = с • ехр(-Т(р - р) / 2),
где с — множитель нормировки (с > о), ¥ — положительно определенная квадратичная форма с матрицей (Ув)-1 (¥(•) > о).
В одномерном случае при нормальном распределении, проверяя статистическую гипотезу, мы имеем в качестве границы критической области две точки, расположенные симметрично относительно центра
предполагаемого распределения.
В четырехмерном же случае, если проверяется гипотеза Н0 : Р = Р*, границей критической области будет эллипсоид Wа с центром в точке Р*. Функция плотности указанного выше распределения принимает на нем постоянное значение, а вероятность того, что точка р с координатами (Ри, Р12, Р21, Р22) принадлежит области, ограниченной этим эллипсоидом, равна 1 -а (а - уровень значимости).
В нашем случае гипотеза Н0 имеет вид Н0 : Р г О. Эту гипотезу можно рассматривать как совокупность простых гипотез Н0 : Р = Р*, где Р* г О. Границей пересечения всех критических областей будет огибающая семейства эллипсоидов Wа, лежащая в О, при условии, что центры эллипсоидов (параметры сдвига семейства распределений) пробегают границу дО. области О. Проверку такой гипотезы можно осуществить по следующему алгоритму, основанному на том, что аргумент и параметр сдвига в функции плотности распределения занимают симметричное положение.
Для распределения Р ~ N(Р, РР) находим такой эллипсоид W с центром в точке Р, чтобы он касался границы дО. области О, а постоянное значение функции плотности этого распределения на границе эллипсоида было бы максимальным. Тогда если вероятностная мера области, ограниченной этим эллипсоидом, будет больше, чем (1 — а), т. е. (1 -у) > (1 -а), или у<а (у - вероятностная мера внешней для эллипсоида части четырехмерного пространства), то гипотеза Н0 может быть отвергнута, и мы, в свою очередь, можем быть уверены в том, что процесс устойчив.
Геометрически задача нахождения эллипсоида W сводится к нахождению точ-
п* 1
ки Р* касания этого эллипсоида с границей дО. области О и имеет стандартное решение по методу множителей Лагранжа [10]. Введем список П := {р, В, V*, Р12, Р21}.
Каждый кусок границы дС1 области О определяется своим уравнением вида ю(р) = 0, где шеП. Каждому выражению ю из списка П соответствует своя функция Лагранжа
4(Р; А) := 1п/(Р) -МР)(шеП),
где ХР) - функция плотности, а уравнением ю(Р) = 0 задается соответствующий кусок границы дО. области О. Из предыдущего раздела видно, что точка Р* есть решение одной из пяти систем уравнений (ш е П) :
(Р; а)
ЗР*
= 0 (/, ] = 1,2), ш(Р) = 0.
За окончательное решение Р* принимаем такое, которое дает большее значение правдоподобия /(Р ) = /, и тем самым получаем искомый эллипсоид
W = {Р : ЯР) = /у},
а вероятностную меру области, ограниченной этим эллипсоидом, сравниваем с выбранным уровнем доверия (1 — а).
В случае относительно небольшого объема выборки для проверки Н0 здесь можно воспользоваться, например, статистикой
г = (Р* -Р)Г^Р-1(Р* -Р),
которая имеет распределение, близкое к Х-квадратному, с числом степеней свободы, равным 4, так как вырожденность РР в этом случае означала бы линейную связь между х и у, т. е. нецелесообразность рассматриваемой модели. Компоненты вектора
я = Рр-1/2(Р* -Р)
здесь будут иметь распределение Стьюден-та с числом степеней свободы V = т — 3.
Кроме х-квадрат-статистики и статистики Стьюдента, для проверки исходной гипотезы можно построить ^-статистику.
Пусть имеет место условие ст1 = ст2, статистическая проверка которого сводится к ^-критерию и является хорошо известным хрестоматийным примером. Тогда исходная гипотеза может быть проверена по-
средством статистики ц ~ г4 2т-6, следует выражению
(р* -р)темет (р* - Р) / 4
которая
Ц :=
М :=
(5! + 52) / (2т - 6)
^ХХТ о ^ о ХХТ
:=! (в,)2(] = 1, 2).
Здесь — суммы квадратов регрессионных остатков уравнений (1о).
Если же условие ст2 = ст2 не проходит статистическую проверку, то для . = 1, 2 имеем статистики ц. как
Р2т-3, вычисляемые
Ц :=
(в* -р.)ТХХТ(в - -в.) / 2 ВД + Р^)/(т - 3) 1
где в* = (в.1, в* 2 )Т, в. = (в. 1, в. 2 )Т.
Независимо от выбора статистики (критерия) вычисление точки касания в* является необходимым, так как перечисленные выше распределения (х-квадрат, Стьюдента и ^-распределение) строятся на основе нормального распределения, а геометрические свойства последнего в этой ситуации играют ключевую роль. Именно эти геометрические свойства и позволяют сводить проверку гипотезы Но :в € О к проверке гипотезы Но : в = в*. Перечисленные выше статистики могут быть ис-
пользованы для проверки любых предположений относительно истинных значений вектора р, т. е. для тестирования модели в целом, а результат тестирования, как было отмечено в самом начале, определяет, согласно общим принципам, правильность выбора модели.
Заключение
Важно отметить, что рассмотренная модель может иметь широкое применение в анализе и прогнозе экономических, социальных и биологических процессов, а также каких-либо иных, которые исследуются в различных областях науки. Все этапы алгоритма построения и анализа этой модели принципиально осуществимы, что, в свою очередь, позволяет рассматривать этот алгоритм в качестве основы для разработки программных продуктов для решения практических задач в различных сферах деятельности.
1. Колемаев В.А. Математическая экономика: Учебник для вузов. М.: ЮНИТИ-ДАНА, 2оо5. 399 с.
2. Прасолов А.В. Динамические модели с запаздыванием и их приложения в экономике и инженерии: Учебное пособие. СПб.: Изд-во «Лань», 2о1о. 192 с.
3. Братусь А.С., Новожилов А.С., Платонов А.П. Динамические системы и модели биологии. М.: Физматлит, 2о1о. 4оо с.
4. Лебедев В.В. Математическое моделирование социально-экономических процессов. М.: «Изограф», 1997. 222 с.
5. Кондрашков А.В. О средствах анализа экономических процессов с использованием непрерывной шкалы времени // Сб. научных статей Всероссийской научно-практической конференции «Проблемы предпринимательства в услови-
СПИСОК ЛИТЕРАТУРЫ
ях кризиса» 3о сентября 2о1о. Т. 1. СПб.: Изд-во СПбГУСЭ, 2о1о. С. 33о-332.
6. Понтрягин Л.С. Обыкновенные дифференциальные уравнения. М.: Наука, 1974. 332 с.
7. Баутин Н.Н., Леонтович Е.А. Методы и приемы качественного исследования динамических систем на плоскости. М.: Наука. Гл. ред. физ.-мат. лит-ры, 199о. 488 с. (Справочная математическая библиотека).
8. Рао С.Р. Линейные статистические методы и их применения. М.: Наука, 1968. 548 с.
9. Себер Дж. Линейный регрессионный анализ. М.: Мир, 198о. 456 с.
10. Пичугин Ю.А. Вероятностные оценки риска и надежности в экономике // Научная сессия ГУАП: Сб. докл.: В 3-х чч. Ч. 1. Технические науки. СПб.: ГУАП СПб., 2о13. С. 195-196.
СВЕДЕНИЯ ОБ АВТОРАХ
КОНДРАШКОВ Анатолий Викторович — кандидат физико-математических наук, доцент кафедры «Прикладные информационные технологии» Санкт-Петербургского государственного экономического университета.
19Ю23, Россия, г. Санкт-Петербург, Садовая ул., 21 ankon69@mail.ru
ПИЧУГИН Юрий Александрович — доктор физико-математических наук, профессор кафедры геометрии Российского государственного педагогического университета им. А.И. Герцена. 191186, Россия, г. Санкт-Петербург, Казанская ул., 6 ушу-рюЬи§т@таП.ги
KondrashkovA.V, Pichugin,Yu.A. ON THE IDENTIFICATION AND STATISTICAL TESTING STABILITY OF THE VOLTERRA MODEL.
Some aspects of the Volterra model are considered in the paper. The emphasis is both on the determination of solution stability area in the space of model parameters and on statistical testing of stability.
volterra model, identification, stability, statistical testing.
REFERENCES
1. Kolemaev V.A. Matematicheskaya ekonomika: Uchebnik dlya vuzov. Moscow, YuNITI-DANA, 2005, 399 p. (rus)
2. Prasolov A.V. Dinamicheskie modeli s za-pazdyvaniem i ikh prilozheniya v ekonomike i in-zhenerii: Uchebnoe posobie. St. Petersburg, Izd-vo «Lan'», 2010, 192 p. (rus)
3. Bratus' A.S., Novozhilov A.S., Platonov A.P. Dinamicheskie sistemy i modeli biologii. Moscow, Fizmatlit, 2010, 400 p. (rus)
4. Lebedev V.V. Matematicheskoe modeliro-vanie sotsial'no-ekonomicheskikh protsessov. Moscow, Izograf, 1997, 222 p. (rus)
5. Kondrashkov A.V. O sredstvakh analiza ekonomicheskikh protsessov s ispol'zovaniem nepreryvnoy shkaly vremeni. Sb. nauchnykh statey Vserossiyskoy nauchno-prakticheskoy konferentsii «Problemy predprinimatel'stva v usloviyakh kriz-
isa» 30 sentyabrya 2010. T. 1. St. Petersburg, Izd-vo SPbGUSE, 2010, pp. 330-332. (rus)
6. PontryaginL.S. Obyknovennye differentsial' nye uravneniya. Moscow, Nauka, 1974, 332 p. (rus)
7. Bautin N.N., Leontovich E.A. Metody i pri-emy kachestvennogo issledovaniya dinamicheskikh sistem na ploskosti. Moscow, Nauka. Gl. red. fiz.-mat. lit-ry, 1990, 488 p. (Spravochnaya matematicheskaya biblioteka). (rus)
8. Rao S.R. Lineynye statisticheskie metody i ikh primeneniya. Moscow, Nauka, 1968, 548 p. (rus)
9. Seber Dzh. Lineynyy regressionnyy analiz. Moscow, Mir, 1980, 456 p. (rus)
10. Pichugin Yu.A. Veroyatnostnye otsenki riska i nadezhnosti v ekonomike. Nauchnaya ses-siya GUAP: Sb. dokl.: V 3-kh chch. Ch. 1. Tekh-nicheskie nauki. St. Petersburg, GUAP SPb, 2013, pp. 195-196. (rus)
THE AUTHORS
KONDRASHKOV Anatoliy V.
St. Petersburg State University of Economics, 21 Sadovaya St., St. Petersburg, 191023, Russia ankon69@mail.ru
PICHUGIN Yuriy A.
Herzen State Pedagogical University of Russia, 6 Kazanskaya St., 191186, St. Petersburg, Russia yury-pichugin@mail.ru
© Санкт-Петербургский государственный политехнический университет, 2014
