Файл: Элементы математического моделирования в программных средах MATLAB 5 и Scilab (Андриевский Фрадков).pdf

ВУЗ: Не указан

Категория: Не указан

Дисциплина: Не указана

Добавлен: 05.04.2024

Просмотров: 432

Скачиваний: 1

ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.

Никто никогда ничего не знает наверняка.

Иосиф Бродский

ГЛАВА 6. П Р И М Е Р Ы РЕШЕНИЯ З А Д А Ч МАТЕМАТИЧЕСКОГО М О Д Е Л И Р О В А Н И Я

6.1.Предсказание курса акций

Вернемся к рассмотренному в п. 3.6 примеру анализа курса акций компании. Рассмотрим задачу предсказания курса по известным данным. В качестве исходного процесса, подлежащего прогнозированию, примем представленный на рис. 3.15, (стр. 101) процесс изменения курса акций компании "Coca Со1ая" после выделения из него линейного тренда [49]. Процесс z[k] имеет нулевое среднее и отличается от исходного процесса у[к] линейной составляющей у[к] = ос + (Зк с параметрами а = 47.3, /3 = 0.0541. Полученное в результате прогноза z[k] используется для предсказания курса акций у[к] по формуле у[к] = у[к] + z[k].

Используем линейную регрессию.

Для рассматриваемых задач наиболее близкой среди из-

вестных линейных моделей

[58] является

авторегрессионная

(АР)-модель. Рассмотрим ее более подробно.

АР-модель представляет собой линейное разностное урав-

нение

 

 

у[к] + аху[к - 1] + а2у[к -

2] + • • • + аПау[к - па] = е[к], (6.1)

где к = 0,1,2,... - дискретное время (номер итерации или

выборочного момента

времени);

а,- (г

= 1,2,...па )

- параме-

тры модели. В это уравнение белый

шум е[к] входит как его

непосредственная ошибка.

 

 

 

Настраиваемые параметры модели образуют вектор

в =

ь а2 , ...,аП в ] .

(6.2)

Введя многочлен от z"1 как A{z)

= l+a!Z"l +a2 z"2 H

banez~n°>

получим описание модели в операторной форме:

 

 

i4(z)y[fc] =

е[к).

 

(6.3)

Приведем уравнения предсказателя для этой модели.

183


Используя (6.1), получаем:

 

 

y[fc|0] = (l-A(z))y[fc].

(6.4)

Теперь введем вектор

 

 

¥>[*] = [ - У [ к ~ 1 ) , .

па]]Т.

(6.5)

Тогда (6.4) можно переписать в виде

 

у[к\6) = <р[к]Т6 =

вГ<р[к],

(6.6)

где вектор параметров модели в определен выражением (6.2).

Таким образом, для систем с одним выходом

предсказатель

описывается в виде скалярного произведения

вектора в и из-

вестного по результатам измерений до момента к - 1 (вклю-

чительно)

вектора <р[к].

 

 

 

В математической статистике модель (6.6)

называют

ли-

нейной регрессией, а вектор

<р - регрессионным

вектором

(ре-

грессором).

Если некоторые из коэффициентов многочлена А

известны

(т.е. не подлежат идентификации), нетрудно полу-

чить линейную регрессию

вида

 

 

 

№ е] = г[к]Те + ф],

 

(6.7)

где слагаемое /х[к] также известно по измерениям до момента

к- 1.

 

 

 

 

При определении параметров моделей нашел широкое при-

менение метод наименьших

квадратов. Лля его применения

уравнения предсказателя (6.7) записываются в виде линейной

регрессии

(6.7).

Ошибка предсказания

6[к}в] =

у[к] - у[к\ в ]

имеет вид 6[к}в]

=

у[к] — (р[к] в. Используем

критериальную

функцию

 

 

 

 

 

 

 

 

JN(6,ZN)

=

-

 

 

<6-8)

 

 

 

 

к=1

 

 

 

Значение

=

argт\перм Jjv(0, ZN)

для данной

функции

Jh(9,Zn) может

быть

получено аналитически;

оно

предста-

вляет собой оценку параметров методом

наименьших

квадра-

тов (МНК, LSE), которая получается из выражения

 

N

(6.9)

к=\

184


Любая оценка, удовлетворяющая уравнениям (6.9), доставляет глобальный минимум функции «7^(0, ZN). Система (6.9) из d линейных уравнений относительно в^ называется системой нормальных уравнений. Если матрица в левой части (6.9) невырожденная, то получается (единственное) значение МНК-оценки Off:

6LNS = I

N

(6.10)

 

 

Уравнения МНК (6.8)—(6.10) можно записать в матричной форме. Лля этого введем УУх^-матрицу

т

 

Г „[!]'

 

т

 

Фл/ =

(6.11)

im1

 

и iVxl-вектор

 

r»w

 

» и

(6.12)

 

ly[N)\

 

Лля многомерного случая, когда у[к]€Ир}

полагаем

^ = [г/[1]Т,1/[2]Т)...,у[ЛГ]т]Т€7г^)

Флг =

Ml],v(2],. ••

~ {Hp* 1)-матрица.

 

Перепишем критерий (6.8)

в виде

 

МО) =

- ф*0||2

= -jfiYN - <S>n0)t(Yn - Ф„в).

(6.13)

Нормальные уравнения (6.9) тогда принимают вид

 

 

 

=

(6.14)

Если det$^<$/v ф 0, получаем МНК-оценку

 

 

eLNs = ( Ф > N)~l$rNYN.

(6.15)

185


Матрица (Ф^Фуу^Фдг является псевдообратной (для Фдг) 1 матрицей, Ф^ = (ФууФуу)~1Флг. Поэтому (6.15) можно переписать в виде

= Ф+ Y„.

(6.16)

Таким образом, уравнение (6.15) дает псевдообратное решение для переопределенной (при N > d) системы линейных уравнений Y^ =

В [58] рассматривается также "взвешенный"

критерий

 

N

2

 

MOtZN)

= £ > [ * ] ( у [ * ] -

 

(6.17)

 

к=\

 

 

или

 

 

 

 

N

2

 

j ^ e ^ » )

= X > W * ] ( y [ f c ] - ф } Т е )

.

(6.18)

к=1

Это связано со следующими причинами:

1.Наблюдения могут иметь различную "надежность". Например, некоторые из них могут содержать большие возмущения и, следовательно, должны иметь меньшие веса.

2.Наблюдения могут иметь изменяющуюся информативность. Например, линейная модель может оказаться слишком грубым приближением для описания системы в некоторой области изменения <р. Наблюдения, относящиеся к такой области, должны иметь меньший вес.

Введем диагональную матрицу Q s = diag{a!,... ,an }. Перепишем критерий (6.17) в виде

 

JN(0) = (Yn - ФNefQN(YN

-

Ф„в) й ||У„ -

 

(6.19)

Минимум (6.19) достигается

при

 

 

 

 

 

On =

 

 

 

 

 

(6.20)

1

Псевдообратной

по

Муру-Пенроузу

матрицей

для всякой

пхга-

матрицы Я называется

матрица

Я+

= \im6 ^o(H

Н + б 2 ! ) ' 1 Я .

Если

Я -

nxn-матрица и det Я ф 0, то

Я +

= Н"1 . Если столбцы Я

линей-

но независимы, то Я +

= (Я Я)~*Я . Выполнено также: (Я ) +

= ( Я + ) ,

Я+ = ( Я Т Я ) + Я Т = Я Т ( Я Т Я ) + , Я Я + Я = Я, Я + Я Я + = Я + (см.,

например,

М ) .

186


Как показано в [58], формулой (6.20) можно пользоваться и для произвольной Qn = Qjsf > 0.

Для критерия (6.18) можно также записать выражение МНКоценки в виде, полностью аналогичном (6.10)

( ж X > w • (6-21)

Будем искать параметры предсказателя для АР-модели вида (6.1)

z[k] + axz[k - 1] + a2z[k - 2] + • • • + anaz[k - n j = е[к]. (6.22)

Уравнения предсказателя имеют вид (6.4), (6.6) с вектором параметров модели в вида (6.2), т.е.

¥>[*]= [-

# - ! ] , . . . , - # -па)}]\

(6.23)

г[к\

9] = <р[к]Те = 6Т<р[к]}

(6.24)

где вектор параметров модели в определен выражением (6.2). Указывается значение р - "сдвиг вперед". Оно показывает, на сколько тактов вперед нужно получить прогноз и значение 71 - порядок АР-модели. По выбранному методу получаются оценки значений параметров АР-процесса. Для этого используется процедура аг, которая входит в состав библиотеки программных модулей SYSTEM IDENTIFICATION ("Идентификация систем") [58, 139, 140]. Возможно использование

методов:

ЧЪ'

- метод прямой и обратной аппроксимации

 

 

(принят по умолчанию);

'Is'

-

метод наименьших квадратов;

'yw'

- метод Юла-Уолкера (Yule-Walker);

'burg'

-

метод Бурга (Burg);

'gl'

-

метод решетчатого фильтра.

Найденные значения параметров регрессионной модели используются при переходе от mema-формы к представлению системы уравнениями состояния: [A,B,C,D,K,X0]=th2ss(th).

187

Здесь параметры А,С,К содержат матрицы регрессионной модели, представленной разностными уравнениями

 

х[к + 1] = Ах[к] + Ке[к],

у[к] = Сх[к]\

(6.25)

матрицы 5, D для АР-модели не используются. Параметр ХО

задает оценку начального состояния системы.

 

В

программе формируется

массив

tm

моментов времени

tk =

кТ0} к = 1,2,... р. Вычисляются

р-я

и (п — 1)-я степе-

ни матрицы А : Ар=А~р; Ап=А~(п-1).

Лалее

для вычисле-

ния

начального состояния

находится

матрица

наблюдаемо-

сти [3, 6, 95, 75]

Г

С?

'

 

 

 

 

Q .

 

 

 

с помощью процедуры obsv тулбокса CONTROL SYSTEMS [74, 61]: Q= obsv (А , С). Эта матрица позволяет найти начальное состояние х0, чтобы получить заданную начальную по-

следовательность выхода {z[-n -hi],... ,z[0]} изт соотношения

xq = Q~lZk> где вектор Zk = [*[fc-n + l], ... ,z[k]] [6]. Процесс формирования вектора Z* и пересчета начального состояния производится на каждом шаге к (начинал с к = п + р) для вычисления прогноза: yl=z(k - p - n+l:k - p); x_=inv(Q)*yl. Из уравнения (6.25) следует, что если известно состояние на А:-м шаге, то его прогноз на fc + p-шаге (полагаем е[к] = 0) получается по формуле х[к + р\к] = Арх[к]. Для вычисления соответствующего значения выхода используем формулу у[к + р\к] =

= САпх[к + р\к] : хр=Ар*х_;

yp(k)=C*An*xp. Лалее

вычисля-

ем ошибку прогнозирования

e = z - y p ( l : l e n g t h ( z ) )

и относи-

тельную среднеквадратическую ошибку

 

 

s e = s t d ( e ) / s t d ( z ) .

Лля восстановления вида исходной последовательности и сравнения ее с прогнозом учитываем линейный тренд

ypr=alpha#ones (size (ур)) +beta*tp' +ур;

Результаты использования АР-модели для заданной последовательности при прогнозе на один и четыре шага представлены на рис. 6.1, 6.2.

Из анализа полученных результатов видно, что точность АР-модели в незначительной степени зависит от выбора метода, а в большей степени - от порядка модели п. С ростом

188