Модель Гудвина: эмпирический анализ в Julia

Julia как инструмент эмпирического анализа моделей Гудвина и Самуэльсона¶

Вводные соображения¶

Предворяя свой эмпирический анализ классических моделей эномических циклов, приведу дословно соображения авторов из статьи "Simulation of business and financial cycles: self-oscillation and synchronization"[1].

"Основной интерес при изучении большинства процессов и явлений, рассматриваемых современной наукой, связан с динамическими закономерностями. Динамика современных технических, биологических, социо-экономических систем и процессов обладает исключительной сложностью.Однако, несмотря на сложность, поведение этих систем на достаточно больших отрезках времени определяют сравнительно простые динамические модели и закономерности, такие, в частности, как предельные циклы, торы, странные хаотические аттракторы, метастабильные со- стояния. ....... Изучению колебательных экономических и социальных явлений посвящено огромное число публикаций, однако, надо признать, что среди них сравнительно мало публикаций, связанных с разработкой и анализом адекватных математических моделей нелинейной динамики социо-экономических систем. Это можно объяснить рядом причин, в частности, сложностью формального описания социально-экономических явлений, трудностью, а, иногда, невозможностью измерения параметров моделей, трудностями экспериментальной проверки моделей, обуслов- ленными экстремально большими периодами колебаний и др."

Если с первым тезисом авторов я могу полностью согласится, то второй требует привести собственные соображения, которые и вызвали к жизни эту работу. Моделей, связанных с математическим описанием экономических циклов, известно достаточно много и многие из них появились довольно давно. При этом теоретический анализ для получения практически значимых выводов часто затруднен по целому ряду причин. Появление в последнее время мощного программно-технического инструментария открывает большие возможности для эмпирического исследования и модификации классических динамических моделей экономических циклов.Большая наглядность, возможность варьирования значительного числа параметров, стохастическая иммитация позволяют эмпирически получать важные теоретические и практические результаты.

В данный небольшой работе я хотел показать на известных моделех Гудвина и Самуэльсона, как в принципе можно осуществлять их модификацию и эмпирический анализ, получая при этом очень наглядные результаты.

Модель Гудвина.¶

Экономические основания модели и ее базовое математическое представление [2]¶

В основе модели лежат переменные, описывающие производство — Y и потребление — C, связанные между собой соотношением:

                 C=αY +β                     (1)

Константа α характеризует пропорциональность между производством и потреблением (α < 1), а константа β — потребление в отсутствие производства (охота, собирательство и т. д.). Еще одна переменная K описывает основной капитал, который включает заводы, оборудование и т. д. Капитал пополняется в процессе производства, что может быть описано уравнением:

                      Y = C + I             (2)

Здесь переменная I, называемая инвестициями, является производной по времени от величины K (далее мы используем стандартное обозначение для производной по времени) :

                  dK/dt = I.                (3)

Исключая переменную C из уравнений (1) и (2), приходим к уравнению:

                  Y =I +β/1−α               (4)

Два уравнения (3) и (4) связывают три переменные: K,Y и I. Третье уравнение может быть получено тем или иным способом из основной идеи модели Гудвина о том, что владельцы капитала через посредство изменения инвестиций стремятся управлять капиталом так, чтобы капитал был адекватен производству. Эта адекватность на математическом языке описывается уравнением:

                  K = γY ,                  (5)

где γ — некоторая константа. В общем случае это означает, что инвестиции следует уменьшать при K > Y и увеличивать в случае противоположного неравенства. Кроме того, следует учитывать определенную ограниченность в стремлении управлять капиталом. Даже если ничего не вкладывать в добавление капитала, он уменьшается вследствие амортизации оборудования. Это означает, что величина I ограничена снизу некоторым отрицательным значением Imin. Из общих соображений можно понять, что эта величина также ограничена сверху некоторым значением Imax. Всё это можно написать в виде математического соотношения:

                   I = f (K −γY ),            (6)

где f — некоторая функция, которая удовлетворяет вышеописанным свойствам, то есть уменьшается с увеличением аргумента и ограничена снизу и сверху. Разные виды моделей, возникающие при том или ином способе моделирования функции f , описаны в книге [7]. Для аналитического исследования уравнений, например, удобно моделировать эту функцию в виде скачкообразного изменения в диапазоне от Imin до Imax, так чтобы между скачками значение I было постоянным. Для численных расчетов, когда решается система дифференциальных уравнений, удобны модели с гладкими функциями, конкретный вид которых будет приведен далее. Из уравнений (3), (4) и (6) следует, что существуют некоторые стационарные значения K0 и Y0, при которых с течением времени все переменные остаются постоянными, причем I0 = 0. Несложные вычисления дают значения:

            Y0 =β/1−α , K0 =γβ/1−α             (7)

Однако, как показывает анализ модели, стационарное состояние является неустойчивым. Если начальное состояние системы характеризуется некоторой точкой вблизи точки стационарного состояния, то с течением времени система отдаляется от стационарного состояния, и начинаются колебания.

Разные виды моделей, возникающие при том или ином способе моделирования функции f , описаны в книге [ ] . Для аналитического исследования уравнений, например, удобно моделировать эту функцию в виде скачкообразного изменения в диапазоне от Imin до Imax так, чтобы между скачками значение I было постоянным. То есть функция f (K −γY ) выбирается в виде:

$$ f\left( {K - \gamma Y} \right) = \left\{ {\begin{array}{*{20}c} {{\mathop{\rm I}\nolimits} _{\min } ,K > \gamma Y} \\ {I_0 ,K = \gamma Y} \\ {{\mathop{\rm I}\nolimits} _{\max } ,K < \gamma Y} \\ \end{array}} \right. $$ Для численных расчетов, когда решается система дифференциальных уравнений удобны модели с гладкими функциями f (K − γY ), конкретный вид которых будет приведен далее. Перед проведением численных расчетов следует еще несколько усложнить модель. Из уравнения (4) следует, что изменения инвестиций мгновенно приводят к изменению производства. В реальности изменения производства несколько запаздывают по сравнению с изменениями инвестиций (см. [7]), так что, вместо соотношения (4), получим:

             Y (t +∆t) = I(t)+β/1−α ,

где ∆t — некоторый малый параметр. Предполагая, что ∆t — малая величина, можно разложить левую часть равенства по степеням ∆t и ограничиться первым членом разложения. В результате получаем уравнение:

        dY(t)/dt ∆t +Y (t) = I(t)+β/1−α     (8)

Аналогично из формулы (6) следует, что изменения величин K и Y мгновенно сказываются на изменении величины I. Однако владельцу капитала необходимо некоторое время ∆t1 для принятия решения об изменении инвестиций. Повторяя аналогичные действия, переходим от уравнения (6) к уравнению:

     dI(t)/dt ∆t1 + I(t) = f(K(t)−γY (t))  , (9)

где ∆t1 — малый параметр. Уравнения (3), (8) и (9) при заданной каким-либо соотношением функции f образуют систему дифференциальных уравнений 1-го порядка, которую можно записать в традиционном виде:

$$ \left\{ {\begin{array}{*{20}c} {\frac{{dK}}{{dt}} = I} \\ {\frac{{dY}}{{dt}} = \frac{1}{{\Delta t}} \cdot \left( {\frac{{I + \beta }}{{1 - \alpha }} - Y} \right)} \\ {\frac{{dI}}{{dt}} = \frac{1}{{\Delta t_1 }} \cdot \left( {f\left( {K - \gamma \cdot Y} \right) - I} \right)} \\ \end{array}} \right. \;\;\;\;(10) $$

Численное исследование этой системы уравнений затруднено слишком большим числом параметров, которые так или иначе нужно задавать при численном моделировании. Уменьшение числа параметров достигается сдвигом и масштабированием переменных.

После целого ряда преобразований можно получить модель Гудвина в виде следующей системы ОДУ:

$$ \left\{ {\begin{array}{*{20}c} {\frac{{dK}}{{dt}} = I} \\ {\frac{{dY}}{{dt}} = \frac{1}{\tau } \cdot \left( {I - Y} \right)} \\ {\frac{{dI}}{{dt}} = \frac{1}{\theta } \cdot \left( {f\left( {K - Y} \right) - I} \right)} \\ \end{array}} \right. \;\;\;\;(11) $$

Здесь $\tau$ и $\theta$ параметры запаздывания влияния I и K. Для дальнейших расчетов нужно выбрать модель для функции f (K −Y ). Она должна убывать с ростом (K − Y) и обращаться в ноль при нулевом значении этой величины. Кроме того, эта функция должна быть ограничена снизу отрицательным значением −µ/(1+µ) и сверху положительным значением 1/(1 + µ). В качестве f можно использовать функция, за основу которой взята функция arctg(x). Подобная функция содержит два параметра, в качестве которых можно выбрать вышеопределенный параметр µ и параметр λ — модуль тангенса наклона кривой в начале координат. При такой параметризации функция f (K −Y ) принимает вид:

$$f\left( {K - Y} \right) = \frac{1}{\pi }\left( {arctg\left( \rho \right) - arctg\left( {\lambda \cdot \pi \left( {1 + \rho ^2 } \right) \cdot \left( {K - Y} \right) + \rho } \right)} \right)\;\;\;\;(12)$$

где $\rho = tg\left( {\frac{\pi }{2} \cdot \frac{{1 - \mu }}{{1 + \mu }}} \right)$.

Численное исследование модели (11)¶

Четыре параметра модели τ,θ,λ и µ имеет довольно ясный смысл. Параметр τ определяет сдвиг по времени между изменением инвестиций и изменением производства.Параметр θ — сдвиг по времени перед принятием решения об изменении инвестиций после анализа состояния системы. Параметр λ — стремление управлять капиталом. Параметр µ — отношение возможной убыли капитала вследствие амортизации оборудования к восстановлению капитала вследствие максимально возможных инвестиций. В соответствии с идеей модели Гудвина, это отношение меньше единицы [7]. Небольшое число параметров позволяет провести некоторое исследование модели на основе численного эксперимента.

В подобных моделях, когда отсутствуют какие-либо соображения по приблизительным величинам параметров и переменных, целесообразно выделить малые параметры, которые положить принимающими значения от 0.1 до 0.9, а остальные параметры положить равными нескольким единицам, или, если параметр по смыслу меньше единицы, — нескольким десятым. Примем простейший набор параметров τ = θ = 0.7, λ = 3.0, µ = 0.35 и при значениях начальных условий в близи нулевой стационарной точки: K0=Y0=I0=0.01. Разработанный Julia скрипт показал, что, в соотвествии с теорией,исследуемые переменные быстро ререходят в колебательнай режим и их атрактор переходит из начальной точки замкнутый ценр.

Углубленный анализ модели Гудвина.¶

Рассмотрим модель экономической динамики Гудвина как модель нелинейного акселератора-мультипликатора, которая может быть записана в виде [3]:

$$\tau \cdot \theta \cdot \frac{{\partial ^2 x}}{{\partial t^2 }} + (\tau + \left( {1 - a} \right) \cdot \theta ) \cdot \frac{{dx}}{{dt}} - \varphi \left( {\frac{{dx}}{{dt}}} \right) + \left( {1 - a} \right) \cdot x\left( t \right) = Q\left( t \right)$$

где $x$ означает доход, $a$ предельный уровень потребления, $τ$ константа представляющая запаздывание в динамике коэффициента развития, $θ$ запаздывание между решением об инвестировании и соответствующими расходами, инвестиции индуцированные изменениями дохода и Q(t) размер независимых издержек по времени. Обобщенная форма этого уравнения была дана Лоренцем:

$$\frac{{\partial ^2 x}}{{\partial t^2 }} + A\left( {x\left( t \right)} \right) \cdot \frac{{dx}}{{dt}} + B\left( {x\left( t \right)} \right) = Q\left( t \right)$$

здесь $x$ – отклонение дохода от равновесия, $A(x)$ – четная функция с условиями $A(0) < 0$; $A′′(0) > 0, B(x)$ – нечетная функция с условием $B(0) = 0$, $Q(t)$ – функция издержек. В своих исследованиях Лоренца и Нуссэ [4] предложили конкретный вид функций $A(x)$ и $B(x)$ :

$$\frac{{\partial ^2 x}}{{\partial t^2 }} + a \cdot \frac{{dx}}{{dt}} \cdot \frac{{\left( {x^2 - 1} \right)}}{{\left( {x^2 + 1} \right)}} - b \cdot x + c \cdot x^3 = Q(t)$$

Функцию $Q(t)$ можно трактовать как некоторое внешнее воздействие. Рассмотрим три случая: $Q(t) = 0$ ; $Q(t)$ -некоторая детерминированная переодическая функция; $Q(t)$ - некоторая стохастическая переодическая функция.

1. Случай $Q(t) = 0$¶

В этом случае модель Гудвина (дифференциальное уравнение второго порядка) можно привести к системе дифуров 1-го порядка:

$${\frac{{dx}}{{dt}} = y}$$ $$\frac{{dy}}{{dt}} = - a \cdot y \cdot {\textstyle{{\left( {x^2 - 1} \right)} \over {\left( {x^2 + 1} \right)}}} + b \cdot x - c \cdot x^3$$

В этой системе существуют три равновесия $M_0(0,0)$ , $M_1 = \left( {\sqrt {\frac{b}{c}} ,0} \right)$, $M_2 = \left( {-\sqrt {\frac{b}{c}} ,0} \right)$.

Равновесие $M_0$ зависит от всех трех параметров. Исследуем эту зависимость. Исходя из экономического смысла параметров: $a$ определяет предельный уровень потребления при соотвествующих значениях $b$ и $c$ ; $0.0< b <1.0$ - норма потребления; $с$-коэффициент снижения дохода, чем он меньше тем больше доход. Зададим следующие значения:$a=[2.0, 2.8]$ ; $b=[0.25, 0.60]$ ; $c=[0.15, 0.2]$. Проведем соответствующие расчеты с помощью разработанного нами скрипта Julia (см. ниже).

Как видно из графиков, в восьми исследованных точках пространства параметров a,b,c мы нашли две точки (2.0,0.25,0.2) и (2.8,0.25,0.20), в которых система из $M_0$ переходит в колебательный режим (аттрактор седло). В остальных исследуемых точках система переходит в устойчивое состояние (аттрактор устойчивый фокус).Если детально исследовать пространство параметров (задать небольшой шаг по каждому параметру,то можно найти точки биффуркации,где система переходит из одного состояния в другое). Например, если в найденных точках параметров, где система находится в колебательном режиме, мы начнем увеличивать значение параметра $а$,то $при a= 3.1839$ система перейдет в устойчивое состояние (точка биффуркации)

Перейдем к анализу устойчивости модели Гудвина в точке $M_1 = \left( {\sqrt {\frac{b}{c}} ,0} \right)$. Зададим $a=2.0$ и посмотрим как изменяется поведение системы в пространстве параметров $b,c$. При этом, из экономических сображений набор значений параметра представим как вектор $B=[0.5,0.4,0.3,0.2]$,а значения параметра $с$ будем определять эксперементально двигаясь от значения параметра $с = 0.28$ вниз с шагом 0.01 с тем, чтобы обнаружить точку бифуркации, где система переходит из одного состояния в другое. Ниже с помощью разработанных нами скриптов, мы получили для каждого значения параметра $b$ серии из 4-х графиков, демонстриющих переход системы из устойчивого состояния в колебательный режим.Полученные результаты можно представить в виде графика(рис. 1). На графике видны точки бифуркации (желтые квадраты),соединяя которые можно получить линию бифуркации.Она делит пространство значений параметров $b,c$ на области устойчивых и неустойчевых состояний.

2. Случай $Q(t)$ детерминированная переодическая функция¶

Пусть внешнее воздействие Q(t) это переодическая функция вида:

$$Q(t) = A_0 \cdot \cos \left( {fr \cdot t + faz} \right) + A_1$$

где $fr$-частота; $faz$-фаза; $A_0$ - амплитуда; $A_1$ - уровень относительно которого происходят колебания.

1. Зададим такие параметры Q(t),при которых она на соотвествующем диапозоне t будет возрастающей функцией. Обеспечить $dQ(t)/dt > 0$ (возрастание функции) можно задав fr=0.02 и faz=10.0. Обеспечить условие Q(t)>0 (экономически это означает наличие внешнего притока средств в экономику) можно выбрав определенное значения $A_1$ и $A_0$ и соотношение между ними.

Пусть множество значений параметра A1 представляет собой вектор [5.0,10.0,20.0,30.0,40.0,50.0,60.0,70.0] (отношение $A_1/A_0$ изменяется от 0.1 до 1.4). Ниже представлен скрипт, который строит серию графиков для каждого значения $A_1$

Из графиков видно, что угол наклона Q(t) увеличивается от arctg(2/30)=3.81 град. до arctg(25/30)=39.81 град. То есть, скорость роста Q(t) увеличивается на интервале исследования более чем в 10 раз. Функция $x(t)$ сначала стабилизируется на определенном уровне, а затем переходит в режим практически линейного устойчивого роста, при этом переходная фаза к усточивому состоянию становиться все менее выраженной (уменьшается число колебаний и их амплитуда). Грубо по графикам оценим скрость роста $x(t)$. При $A_1/A_0=0.1$ $\Delta x \approx 0.0$ ; $A_1/A_0=0.2$ $\Delta x \approx 0.2$; $A_1/A_0=0.4$ $\Delta x \approx 0.4$; $A_1/A_0=0.6$ $\Delta x \approx 0.8$; $A_1/A_0=0.8$ $\Delta x \approx 1.6$; $A_1/A_0=1.0$ $\Delta x \approx 3.0$; $A_1/A_0=1.2$ $\Delta x \approx 4.6$. В точке, где $A_1/A_0=1.2$ кривая роста стновитья явно выпуклой и это последняя точка, где на всем диапозоне исследования $x(t)>0$ ($x(0)\approx 0.0$). В точке $A_1/A_0=1.4$ $x(t)$ и $Q(t)$ уходят в отрицательную зону и движение $x(t)$ становиться хаотичным. Точкой бифуркации является точка расположенная между $A_1/A_0=1.0$ и $A_1/A_0=1.2$. Это хорошо видно по фазовым портретам последних двух точек исследования.

2. Зададим такие параметры Q(t),при которых она на соотвествующем диапозоне t [0,100] будет убывающей функцией. Обеспечить $dQ(t)/dt < 0$ (убывание функции) можно задав fr=0.02 и faz=0.0. При этом на исследуемом инттервале времени функция может $Q(t)>0$ и $Q(t)<0$ (экономически это означает наличие убывающего притока средств в экономику, переходящего в отток ресурсов из экономики ).Данное поведение $Q(t)$ можно обеспечить, выбрав определенное значения $A_1$ и $A_0$ и соотношение между ними.

Пусть множество значений параметра A0 представляет собой вектор [100.0, 50.0, 40.0,35.0] (отношение $A_1/A_0$ изменяется от 1.0 до 0.35). Ниже представлен скрипт, который строит серию графиков для каждого значения $A_0$

Из графиков видно, что по мере движения от значения отношения $A_1/A_0 =1.0$ до $A_1/A_0 =0.4$ система деградирует (x(t) падает) со все увеличивающейся скоростью, а между точками 0.4 и 0.35 она коллапсирует и разрушается

3. Зададим такие параметры $Q(t)$,при которых она на соотвествующем диапозоне t [0,100] будет переодической функцией ($fr=0.2 , faz=0.0$). При этом на исследуемом инттервале времени функция может $Q(t)>0$ и $Q(t)<0$ (экономически это означает наличие переодически увеливающегося и убывающего притока средств в экономику, переходящего в определенные моменты t в отток ресурсов из экономики).Данное поведение $Q(t)$ можно обеспечить, выбрав определенное значения $A_1$ и $A_0$ и соотношение между ними.

Пусть множество значений параметра A0 представляет собой вектор [200.0, 10.0, 2.0] (отношение $A_1/A_0$ изменяется от 0.025 до 2.5). Ниже представлен скрипт, который строит серию графиков для каждого значения $A_0$

Из графиков видно, что при значении отношения $A_1/A_0 =0.025$ система сохраняет устойчивое состояние(x(t)≈ 9.2 ).При $A_1/A_0=1.0$ находится в колебательном режиме $1.8 ≤ x(t) ≤ 4.2$ с частотой близкой к частоте $Q(t).При $A_1/A_0 =2.5$ система переходит в сложное состояние с переодической составляющей с постоянным переходом в отрицательную зону,с экономической точки зрения это соотвествует коллапсу и разрушению системы (фазовый портрет седло с петлями).

Q(t) стохастическая функция¶

Используем в качестве Q(t) переодическую функцию, в которой вносятся случайеые положительные возмущения в частоту и aмплитуду. Ниже представлен скрипт, который генерирует $N$ реализаций подобной функции

Статистический анализ Q(t)¶

На множестве,состоящем из $N$ реализаций $Q(t)$, в каждом временном сечении сформирована случайная величина, то есть имеется система из T случайных велечин (в нашем случае T=20). В скрипте определяются среднее (mean),среднеквадратичное отклонение (std) и коэффициет вариации (Cv) для каждой случайной величины q(t),t=1:20. На этой основе строятся графики:mean(t),mean(t)+3std(t),mean(t)-3std(t) и Cv(t). Из графиков видно, что в 99% реализаций Q(t) находятся в интервале [0.0,105.0]. Минимум коэффициета вариации min Cv = 18.32 %, а максимум max Cv=32.55. Посмотрем как Q(t) c такими статистическими характеристиеами влияет на x(t) и каков при этом фазовый портрет, порождаемый моделью Гудвина