Рассматривается задача идентификации динамической модели рынка вальрасовского типа со многими товарами, описываемой системой дифференциальных уравнений с запаздывающими аргументами, разрешенных относительно производных, на основе метода наименьших квадратов. Получена вычислительная структура алгоритма идентификации неизвестных коэффициентов системы и запаздываний. Работа алгоритма иллюстрируется модельными численными примерами и примером идентификации реального рынка одного из видов компьютерных комплектующих - видеокарт различных производителей.
THE IDENTIFICATION OF DYNAMIC MARKET MODELS OF WALRASIAN TYPE WITH MANY GOODS .pdf Идентификация моделей рынка вальрасовского типа, описываемых дифференциальными уравнениями с запаздывающим аргументом, осложнена наличием временных лагов, запаздываний реакций поставщиков товаров на изменение цен этих товаров. Работ по идентификации дифференциальных уравнений с запаздывающим аргументом довольно мало. Среди них можно указать, например, работы D.W. Brewer [1], C. Baker и E. I. Parmuzin [2], В.И. Ловчакова [3], А.В. Прасолова [4], В.Ф. Лебедева и Е.А. Ситникова [5], С.А. Минюка и А.В. Метельского [6].В данной работе предлагается конструктивный алгоритм идентификации параметров динамической модели рынка многих товаров (коэффициентов модели и запаздываний) на основе метода наименьших квадратов. Частный случай идентификации рынка одного товара рассмотрен в работе [7].1. Алгоритм идентификации многомерной динамической моделиРассмотрим динамическую модель рынка вальрасовского типа с N товарами, описываемую системой дифференциальных уравнений с запаздывающими аргументами [8]:Т^-Y" ^*(t))"Q(P(tТ))]"-X 5=1, yy[Qf (Pj (t))- QSj (Pj (t- T j))] , * e [t0, T ],P(t) = p(t) , t e [t0-ti,t0), Pi(t0) = Po , i = 17N, (1) где P (t) - рыночная цена i-го товара, xt - задержка (временной лаг) поставки i-го товара на рынок, QD (P (t)) и gf (р (t -xt)) - объемы спроса и предложения i-го товара в момент времени t, р (t) - начальная функция (цена i-го товара до начального момента времени t0), i = 1, N.P , Qi ) рыночных равновесий по Вальрасу определяются пересечением линий спроса и предложения [9]: Q/3 (р) = Q* -at (p -p*), Qi (P) = Qi + ßi (PJ _ P ), i = 1, N, где P - равновесная цена i-го товара, Qi = Qi (Pi ) = Qi (Pi ) - равновесный объем спроса-предложения на i-й товар,* Ii*\ *a > 0 , ßf = Qi (Pf - Pmin i) > 0 , Pmin i < Pi - минимальная цена, по которой продавец соглашается продавать i-й товар. Коэффициенты угу > 0 характеризуютконкуренцию товаров на рынке. Введем обозначения:P(t) = (P(t) P2 (t) ... Pn (t) )T,P{t-T) = (P (t) P2 (t-T2 ) ... Pv (t-Tn ))T ,* ч TP = (P P2 ... Pv ) ,где у =dtdtА = уа ,B = yß ,Taj0...0a2 ...00...d In PN (t) dt(ß 00 ß2000 ^С учетом этих обозначений, подставив gf (р (г)) и gf (р (t -тг-)) в (1), приве-дем систему (1) к видуd ln P(t)dtA(P(t) - P ) + B(P(t -t) - P ),(2)где параметры A, B, P и т, вообще говоря, неизвестны и подлежат определению по результатам наблюдений за ценами товаров на некотором интервале времени to < t < T . *Для идентификации неизвестных параметров A, B, P и т по наблюдаемому решению P (t) уравнения (2) воспользуемся методом наименьших квадратов.Проинтегрируем формально уравнение (2) на интервале [t0, t]:ln I AÄ L A J (P (t) - P* )dt-l(t > t0) +B ((P (t) - P* )dt-т> t0), (3)где 1(.) - индикатор события, указанного в скобках. Разобьем интервал интегрирования [t0, T] на n примыкающих друг к другу интервалов (tk, tk+1 ], k = 0, n -1, длины h = (T -10 )ln .При t = tk равенство (3) примет видПри t = tk+l:ln/£Ы1 = a)(P(t)-P*)dtA(tk >tQ) +V P0 J to+B"J (P(t)-P*)dt-\(tk-т>to).to^p1-1 = A^f (P(t)-P*)dt+1 >toV(4)+B J (P(t)-P*)dt-т>to) . На каждом k-м интервале вычтем равенство (4) из (5):ln2V P l'i ) У tk +Bt+) \p(t)-P*)dt-т>to) . Считая шаг h малым, представим интегралы в (6) по формуле трапеций: P (tk+I )^ hA[ P (tk+i) + P (tk)-2P* ]_P (tk)l(tk ^t0) +(5)(6)2ГP (tk+1 -x) + P (tk -x)-2P* 1(7)Введем следующие обозначения:ln | pN (tk±1 )ln|Äifk±l)l ln ( P (tk+l )^ Г1 Pl ('к ) ) "\ P2 (tk ) ) "' ~l Pn (tk )Фк(P*) = h(P(tk+1)-2P* + P(tk))-1(tk > t0)/2 , Vk(P*,t) = h(P(tk+1 -t)-2P* + P(tk -x)) 1 (tk - x > t0)/2 ,где Xk , Фк (P*), V|/k (P*,T) _ векторы-столбцы размерности N. Тогда равенство (7) примет видXk = (P*) + (P*, x), k = 0, n-1.В матричной форме это соотношение запишется так:x = Y (P*, т) C,матрица коэффициентов размерности(8)(9)(10)2 N х N,Y (P*, т) = (ф (P*) у (P*, т)) - матрица плана размерности n х 2 N , ах, ф, уцатрицы размерности n х N : х = (x0 хх xn-l )T, ф = (ф0 ф1 фп-1 )T ,V = (Vo Vi Vя-1 )T .МНК-оценки неизвестных параметров должны обеспечить минимум суммы квадратов невязок системы равенств (10):J = tr((x-Y(p*,t)c)T (x-Y(p*,t)c)) => min, . (11)Находя отсюда МНК-оценки матриц коэффициентов A и BС:(У (P*, т)T Y(P*, т)) Y(P*, т)T x (12)через P , т и подставляя их в выражение для J, получим целевую функцию для оценивания P*,т:J = tr(x- Y(P*,t)(y(P*,Y(P*,t)) 1 Y(P*,tfx | xvvx| x-Перебирая с некоторым шагом значения P* и т, находим точку P*, т минимума J(P*, т) и вычисляем C . В случае, когда в качестве наблюдений используются решения системы, полученные с помощью метода шагов по схеме Рунге-Кутты 4-го порядка точности, ошибки, с которыми вычисляются параметры модели, - это ошибки квантования, возникшие за счет рассогласования схем второго и четвертого порядков точности.2. Идентификация динамической модели рынка вальрасовского типас одним товаромРассмотрим модель рынка вальрасовского типа при N = 1. Пусть до момента времени t0 рынок находился в состоянии равновесия, а в момент t0 цена товараскачком изменилась от значения P* до значения, например, P0 = 2P*. Возникает переходный процесс, в результате которого рынок либо постепенно снова приходит к состоянию равновесия, если т< ткр - запаздывание меньше критического,либо устанавливаются циклические колебания цены, если т> ткр . Численное решение начальной задачи для нелинейного дифференциального уравнения с запаздывающим аргументом (2) с указанными начальными условиями легко получить методом шагов [10] с использованием вычислительной схемы Рунге-Кутты 4-го порядка точности [11].При t0 = 0, T = 100 и значениях параметров P* = 1, Q* = 1, у = 0,7, а = 0,5 , Pmin = 0,1P*, так что ß = Q'liP* - Pmin) = 10/9 , A = -0,35 , B = -7/9 = -0,77778, решение показано на рис. 1 и 2 при т = 1,5 < ткр (положение равновесия P* устойчиво) и на рис. 3 и 4 при т = 3,5 > ткр (положение равновесия P* неустойчиво). Шаг интегрирования метода Рунге-Кутты будем брать равным h = 10-3.По полученным значениям P(t) идентифицируем с помощью описанного метода параметры A, B, P*, т. Шаг идентификации (шаг отсчета наблюдений) возьмем равным h = 10-3.На рис. 5 - 8 представлен вид поверхности J (P, т) и ее горизонтальные и вертикальные сечения при истинных P* = 1 и т0 = 1,5 < ткр (равновесие устойчиво), а на рис. 9 - 12 - то же, но при т0 = 3,5 > ткр (равновесие неустойчиво). Хорошо видно, что целевая функция МНК J (P, т) достигает минимума по P и т именно при этих значениях параметров. J-103//1,5 -/1 -/0,5 -1101 2 3 4 5 Рис. 10. Линии уровней /(Р,т)при То = 3,5 > ТкрВ табл. 1 приведены значения оценок параметров A и B модели, значения целевой функции и разница между реальными значениями и оценками параметров A и B при т = 1,5 и т = 3,5.Таблица 1 Результаты идентификации при n = 1Параметрыt 1,53,5A-0,349966-0,349966B-0,777768-0,777742j5,61x10-131,21x10""a - Al3,43 x10-53,4 x10-5B - B1,01x10~53,55 x10-5Значение минимума функционала при т = 3,5 больше, чем при т = 1,5, оценки параметров получаются немного хуже, чем при т = 1,5. Можно показать, что погрешность идентификации минимальна при совпадении шага идентификации с шагом интегрирования.3. Идентификация динамической модели рынка вальрасовского типас двумя товарамиРассмотрим модель рынка вальрасовского типа с двумя товарами. Пусть аналогично предыдущему пункту до момента времени t0 рынок находился в состоянии равновесия, а в момент t0 цены обоих товаров резко изменились:P (to) = 2P*, P2 (to) = 1.5P*.Численное решение начальной задачи для нелинейного дифференциального уравнения с запаздывающим аргументом (2) с указанными начальными условиями получим, как и в предыдущем пункте, с использованием метода шагов и схемы Рун-ге-Кутты 4-го порядка точности.Рассмотрим следующие значения параметров:q$$$$$ $-0,6 0,2 ^ (0,9^ р 0,3 -0,5J ' ""11,3J ' Pmß=f10/9l, ^=f-°'54 0,26 1, *=1.10/6j I. 0,27 -0,65Jh = 10-3, p = p = P = 1, Q = g2 = Q = 1,у:'0,1 ^ .0,4 Утак что-6/9 1/3 ^ 1/3 -5/6^' Решив систему, получим некоторое решениеP(t ) =IP2 (t)( P (t Упоказанное на рис. 13 и 14 при 1 = 15! (положение равновесия P устойчиво) ина рис. 15 и 16 при т =(положение равновесия P неустойчиво).0,5Q si (t), Q 2 (t)Рис. 13. p (t), P2(t) и (t), Q2S (t),Рис. 14. QXS (p) и ß/ (P2)зависимость P2 (t) от p (?) при т = (3 5)при т = (3 5)По полученным значениям P(t) идентифицируем параметры модели A, B, т. Шаг идентификации возьмем равным 10 . Значения параметров р = P2 = P = 1будем считать известными.На рис. 17 - 20 представлен вид поверхности J(т) и ее горизонтальные и вертикальные сечения при истинных значениях т1 = 3 и т2 = 5 (равновесие устойчиво), а на рис. 21 - 24 - то же, но при т1 = 6 и т2 = 7 (равновесие неустойчиво). Хорошо видно, что целевая функция МНК J (т) достигает минимума по т1 и т2именно при истинных значениях параметров.Поверхности J(т), полученные для модели рынка вальрасовского типа с двумя товарами, обладают интересной особенностью: рассматривая сечения поверхности при любом фиксированном значении одного из запаздываний, получаем, что кроме минимума, достигаемого в точке второго запаздывания, кривая сечения имеет излом в точке, соответствующей значению, в котором достигается минимум по первому запаздыванию. Например, на рис. 19, где представлены сеченияJ (т) как функции т1 при различных т2 при т = (3 5)T, минимум достигается при т1 = 3 , а излом - при т1 = 5 . На рис. 20, где представлены сечения J(т) какфункции т2 при различных т1 при т = (3 5)T, минимум достигается при т2 = 5 , а излом можно наблюдать при т2 = 3 . Таким образом, на поверхности J (т) можно наблюдать четыре «борозды»: две из них проходят параллельно оси т2, две -параллельно оси т . «Борозды» соответствуют значениям т = 3 , т = 5 и т2 = 3 , т2 = 5 . В одном из пересечений таких «борозд» и находится минимум поверхности. В случае т = (6 7)т «борозды» более сглажены и не так заметны.TВ табл. 2 показаны значения оценок параметров a и b модели, значения целевой функции и разница между реальными значениями и оценками параметров a иb при т = (3 5)T и т = (6 7)T .В таблице под |a - AI, \ß -ß\ понимается следующее:A - A =Ич AniA12 A12|(\B11- B11lB12 - B1'2122A22 AB22 B22При больших значениях запаздываний т = (6 7) модель рынка идентифицируется немного хуже, чем при т = (3 5)T . Можно показать, что в рассматриваемом случае N = 2, как и в случае N = 1, погрешность идентификации будет минимальной при совпадении шага идентификации с шагом интегрирования.4. Идентификация модели вальрасовского типа на примере рынка компьютерных комплектующихРассмотрим теперь решение задачи идентификации динамической модели рынка вальрасовского типа по реальным данным динамики цен на рынке компьютерных комплектующих. Возьмем для примера три товара - видеокарты разных производителей, конкурирующих между собой:1))Видеокарта 128Mb DDR Gigabyte GV-RX55HM256P8-RH (OEM) 64bit+DVI+TV Out .2))Видеокарта 128Mb DDR ASUSTeK EAX550HM512/TD (RTL) 64bit+DVI+TV Out .3)Видеокарта SVGA 128Mb DDR Sapphire (OEM) 64bit +DVI+TV Out.Наблюдения за ценами этих товаров взяты в интернет-магазине электронной техники «Техноград» [12]. Временной интервал наблюдения составляет 132 дня с 22.11.2006 по 03.04.2007 года.Пусть гнабл1, 4абл2, 4аблз - векторы моментов времени, в которые были сделаны наблюдения за ценами соответственно 1-го, 2-го и 3-го товаров; Рнабл1, Рнабл2, Рнабл3 - векторы наблюдаемых значений цен этих товаров; пнабл1, пнабл2, пнабл3 - объемы наблюдений, т.е размерности векторов моментов времени и наблюдаемых значений цен товаров. Динамика наблюдений за ценами рассматриваемых товаров изображена на рис. 25.1600 1500 1400 1300 12001100Рис. 25. Наблюдения за ценами Ртбл1, Рнабл2, РнаблЗ видеокарт разных производителейОчевидно, что равновесные цены рассматриваемых товаров различны и изменяются с течением времени. Для того чтобы можно было воспользоваться для рынка комплектующих моделью вальрасовского типа и идентифицировать ее по данным наблюдений, необходимо исключить тренд равновесных цен. Кроме того, удобно привести наблюдения к безразмерному виду. Тогда равновесные цены всех товаров будут одинаковыми и равными единице. Для этого отдельно для каждого товара по имеющимся наблюдениям выделим линейный тренд видаР* = ajIj + bj/наблj , j = 1,2,3 , (14)где Pj* - инаблу - вектор-столбец значений равновесной цены у-го товара в моменты времени наблюдений, a,, bj - неизвестные параметры тренда, скаляры, подлежащие определению по результатам наблюдений, - единичный вектор-столбец размерности пнабл,.Для нахождения оценок неизвестных параметров aj, bj, j = 1,2,3, воспользуемся методом наименьших квадратов:"наблj 2Jj = X (Рнабл ji - (aj + bj'набл ji)) => mm , j = 1> 2, 3 ,i=1a j,b jили, в векторной форме,j = \pj I набл j2min, j = 1,2,3 .(15)Введем матрицы плана и векторы неизвестных параметров тренда для каждого из товаров:Yj (1 'набл j ) 'Cn.\v bj j, j= 1,2,3 .Тогда (15) перепишется следующим образом:Jj = (Рнабл j - Yj cj f (Рнабл j - Yj cj ) => min , j = 1, 2, 3 .Отсюда находим МНК-оценки параметров:(:TYj У' YjTРнаблj, j = 1,2,3 .На рис. 26 - 28 показаны наблюдаемые значения цен и полученные линейные тренды - равновесные цены соответственно первого, второго и третьего товаров.13001600Рис. 26. Наблюдения за ценой Рнабл1 и равновесная цена PJ* = 1,48 х 1031Х -1,324?на6л113001600йчРис.28. Наблюдения за ценой Рнаблз и равновесная цена P3* = 1,565 х 10313 - 2,467tHПриведем наблюдения к безразмерным единичным ценам следующим обра-зом:пр jiнабл jiji7 = 1,2,3 , i = 0,nH1.Рассмотрим временной промежуток [t0, T] длиной 132 дня, где t0 = 0,T = 132. Выбрав достаточно малый шаг интегрирования h = 10-3, обеспечивающий приемлемую точность решения задачи идентификации, разобьем промежуток [t0 ,T] точками tk+1 = tk + h, k = 0, n -1, на n = (T -10)/h примыкающих друг к другу интервалов одинаковой длины h. На каждом промежутке [tнабл i' набл i+1 ] ,i = 0, пнабл j -1, определим интерполяционные значения наблюдаемых цен товаров в точках tk:jmP- PDнабл ji+1 набл ji , >Рнабл ji(tm ^набл ji ) ,наблнабл jij = 1,2,3 , m = гнаблij h, гнаблi+ljh .По полученным таким образом значениям цен товаров р (t), P2 (t), P3 (t) будемпроводить идентификацию модели вальрасовского типа рынка компьютерных комплектующих.На рис. 29 показаны цены рассматриваемых товаров, приведенные, как описано выше, к единичной равновесной цене.а.0Задача идентификации решалась по формулам (12) - (13) при фиксированном значении P* = 1. При идентификации получены следующие значения параметров модели:(14,56 1,43 11,37 ^f-14,58 -1,39 -11,37^ (0,0Г\A =-0,46 8,991,6 -0,66 -2,39 7,73B =0 47-8,96 -1,61 2,22 0,67-7,71 0,010,01Запаздывания т получились ничтожно малыми. И действительно, на рассматриваемом рынке компьютерных комплектующих время доставки товара на рынок составляет 3 - 7 дней, тогда как цена за это время практически не изменяется, поэтому такой рынок по сути не содержит запаздываний, колебания цен происходят по другим причинам, одной из которых является наличие конкуренции между товарами.На рис. 30 - 32 приведены для сравнения наблюдаемые приведенные цены товаров P2W, Рз(0 и траектории цен Pvan(t), Pvaa(t), Рш/з(0, полученные при решении идентифицированных по параметрам A, B, т дифференциальных уравнений, описывающих динамику рынка. При решении предполагалось, что до момента наблюдений t0 = 0 рынок находился в состоянии равновесия.соа?Возврат к исходным ценам товаров проводится по обратным формулам:Pji = Pvalji Pji ' J => i = 0> пнабл j - 1 Полученные траектории в сравнении с исходными наблюдениями за ценами товаров приведены на рис. 33 - 35.1600® 1500I1300S 14001600 1500яS 1400
Интернет-магазин электронной техники «Техноград». [Электронный ресурс]. Режим доступа: <http://technograd.tomsk.ru/>, свободный.
Эльсгольц Л.Э. Введение в теорию дифференциальных уравнений с отклоняющимся аргументом. М: Наука, 1964. 128 с.
Тихонов А.Н., Васильева А.Б., Свешников А.Г. Дифференциальные уравнения. М.: Наука, 1980. 232 с.
Гальперин В.М., Игнатьев С.М., Моргунов В.И. Микроэкономика: В 2 т. / Общ. ред. В.М. Гальперина. СПб.: Экономическая школа, 2002. Т. 1. 349 с.
Поддубный В.В., Сухарева Е.А. Исследование динамической модели рынка вальрасовского типа со многими товарами // Вестник ТГУ. 2006. № 293. С. 53 - 58.
Сухарева Е.А. Идентификация динамической модели рынка вальрасовского типа // Научное творчество молодежи: Материалы XI Всеросс. науч.-практич. конф. (20 - 21 апреля 2007 г.). 2007. С. 50 - 54.
Минюк С.А., Метельский А.В. О построении непрерывной восстанавливающей операции в задаче полной идентификации линейных стационарных систем с запаздыванием // Дифф. уравн. 2003. Т. 39. № 8. С. 1052 - 1057.
Прасолов А.В. Математические модели динамики в экономике. СПб: Изд-во СПбГУ, 2000. 300 с.
Лебедев В.Ф., Ситников Е.А. Идентификация систем с запаздыванием с использованием полиномиальных рядов // Системы управления и информационные технологии: Межвузовский сборник научных трудов. Вып. 9. Воронеж: Центрально-Черноземное книжное издательство, 2002. С. 61 - 65.
Ловчаков В.И. Идентификация линейных динамических систем с запаздыванием // Дифференциальные уравнения и их приложения. Тула: ТулГУ, 1996. С. 37 - 44.
Brewer D.W. Gradient methods for identification of distributed parameter systems // Decision and Control: Proc. of the 28th IEEE Conf. 1989. V. 1. P. 599 - 603.
Baker C., Parmuzin E.I. Analysis via integral equations of an identification problem for delay differential equations // J. Integr. Equat. and Appl. 2004. V. 16. No. 2. P. 111 - 117.