Имитация реального делового цикла

По сути, мне нужно воспроизвести «Руководство пользователя по решению моделей реального бизнес-цикла» Хартли ( http://www.econ.ucdavis.edu/faculty/kdsalyer/LECTURES/Ecn235a/Linearization/ugfinal.pdf ). В частности, я хочу смоделировать динамическую систему, подразумеваемую моделью, которая указана следующим образом:

где - потребление, - предложение рабочей силы, - капитал, - авторегрессивный технологический процесс, - выпуск, а - инвестиции.$c$$h$$k$$z$$y$$i$

Я моделирую это, используя следующую логику: скажем, в момент времени все находится в устойчивом состоянии, и все значения равны 0, из чего мы получаем . Затем, в момент времени , давая шок системе через , я для и (так как у меня есть «шокированный» и полученный ранее . Затем я подключаю эти два для извлечения остальных, а именно - и повторяю процесс.$t$$k_{t+1}$$t+1$$\varepsilon$$c_{t+1}$$h_{t+1}$$z_{t+1}$$k_{t+1}$$y_{t+1}, i_{t+1}, k_{t+2}$

К сожалению, я получаю взрывной процесс, который не имеет смысла:

Я также включил код R, который используется для симуляции этого:

n<-300

data.simulated <- data.table(t = 0, zval = 0, cval = 0, hval = 0, kval = 0, yval = 0, ival = 0)
data.simulated <- rbind(data.simulated, data.table(t = 1, kval = 0), fill = TRUE)

for (ii in 1:n){

  ##initial shocks
  eps <- rnorm(1, mean = 0, sd = 0.007)
  zt1 <- data.simulated[t == ii - 1, zval]*0.95 + eps
  kt1 <- data.simulated[t == ii, kval]

  ##solve for ct, ht
  lmat <- matrix(c(1, -0.54, 2.78, 1), byrow = T, ncol = 2)
  rmat <- matrix(c(0.02 * kt1 + 0.44 * zt1, kt1 + 2.78 * zt1), ncol = 1)

  solution <- solve(lmat, rmat)
  ct1 <- solution[1, ]
  ht1 <- solution[2, ]

  ##now solve for yt1 and kt2 and it1
  yt1 <- zt1 + 0.36 * kt1 + 0.64 * ht1
  kt2 <- -0.07 * ct1 + 1.01 * kt1 + 0.06 * ht1 + 0.1 * zt1
  it1 <- 3.92 * yt1 - 2.92 * ct1

  ##add to the data.table the results
  data.simulated[t == ii, c("zval", "cval", "hval", "yval", "ival") := list(zt1, ct1, ht1, yt1, it1)]
  data.simulated <- rbind(data.simulated, data.table(t = ii + 1, kval = kt2), fill = TRUE)
}


a <- data.simulated[, list(t, cval, ival, yval)]
a <- data.table:::melt.data.table(a, id.vars = "t")
ggplot(data = a, aes(x = t, y = value, col = variable)) + geom_line()

Мой вопрос прост - система, указанная в статье, по своей сути нестабильна и дает результаты, или я где-то допустил ошибку?

взрывчатость

Статья содержит ошибку, которая вызывает взрывную динамику в вашем моделировании (хотя предположительно основные расчеты в статье были правильными). Условие равновесия, полученное из разложения по собственным значениям, содержится в третьей строке матрицы на странице 12 статьи, переменные упорядочены как (я опущу тильды, так что все строчные переменные следует понимать как логарифмические отклонения). Сравнивая с уравнением (16) на с. 13 мы видим, что коэффициенты для и переключаются, и поэтому правильное условие$Q^{-1}$$(c,k,h,z)$$k$$h$

$$c_t = 0.54 k_t + 0.02 h_t + 0.44 z_t$$

моделирование

Во-первых, мы можем выразить потребление и труд в виде линейной функции переменных состояния (нет необходимости решать систему на каждом этапе моделирования). Условия межвременного и внутривременного равновесия можно записать в виде

$$\begin{bmatrix}1 & -0.02 \\ 2.78 & 1 \end{bmatrix} \begin{bmatrix} c_t \\ h_t\end{bmatrix} = \begin{bmatrix} 0.54 & 0.44 \\ 1 & 2.78 \end{bmatrix} \begin{bmatrix} k_t \\ z_t\end{bmatrix}$$

поэтому после умножения на обратное мы получаем

$$\begin{bmatrix} c_t \\ h_t\end{bmatrix} = \begin{bmatrix} 0.53 & 0.47 \\ -0.47 & 1.47 \end{bmatrix} \begin{bmatrix} k_t \\ z_t\end{bmatrix}$$

Далее переход для состояний можно записать как

$$\begin{bmatrix} k_{t+1} \\ z_{t+1} \end{bmatrix} = \begin{bmatrix} -0.07 & 0.06 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} c_t \\ h_t\end{bmatrix} + \begin{bmatrix} 1.01 & 0.1 \\ 0 & 0.95 \end{bmatrix} \begin{bmatrix} k_t \\ z_t\end{bmatrix} + \begin{bmatrix} 0 \\ \epsilon_{t+1}\end{bmatrix}$$

который может быть уменьшен путем замены контрольных переменных на

$$\begin{bmatrix} k_{t+1} \\ z_{t+1} \end{bmatrix} = \begin{bmatrix} 0.94 & 0.16 \\ 0 & 0.95 \end{bmatrix} \begin{bmatrix} k_t \\ z_t\end{bmatrix} + \begin{bmatrix} 0 \\ \epsilon_{t+1}\end{bmatrix}$$

Теперь симуляция должна быть тривиальной, вот пример Matlab / Octave:

T = 200;
X = zeros(2,T);
for i=2:T
    X(:,i) = [0.94 0.16; 0 0.95] * X(:,i-1) + [0; 0.007*randn()];
end
Y = [0.53 0.47; -0.47 1.47] * X;
figure
plot(1:T, [X; Y])
legend('k','z','c','h')

Конечно, на практике вам, вероятно, следует пересчитать все решение в целом, включая разложение по собственным значениям, чтобы вы могли изменять параметры и т. Д.

Заключительные новости 20 марта 2015 года : Я получил письмо от проф. К. Салер, один из авторов Руководства пользователя. В повторном сообщении он подтвердил, что обе проблемы (см. Мой ответ ниже, а также ответ @ivansml) существуют:

а) Правильное уравнение для закона движения потребления, как показывает @ivansml

б) Число будет ошибочно называют «дисперсией» (стр. 11) в статье. На самом деле, это стандартное отклонение , и, действительно, такая величина является типичной находкой в ​​данных (проф. Сэлиер ссылается на «стр. 22 формы 1 Кули и Прескотта« Исследование границ делового цикла »).$0.007$

Обе ошибки относятся к печатной версии статьи. Другими словами, моделирование за рисунком 1 документа является правильным: они используют правильное уравнение для потребления, и они используют в качестве стандартного отклонения возмущения в технологическом процессе. Так что это второй график ниже, что он действителен.$0.007$

---

ФАЗА A
Я с помощью моделирования (и с использованием правильного стандартного отклонения) проверил, что модель взрывается, хотя она делает это вверх, а не вниз. В статье должна быть расчетливая ошибка, которая, тем не менее, почему-то не была «передана» для моделирования авторов. На данный момент я не могу думать ни о чем другом, так как методология стандартная. Я заинтригован и до сих пор работаю над этим.

ФАЗА B
После ответа @ ivansml, который выявил ошибку в статье (которая, по-видимому, не была сделана при моделировании авторов) , я думаю, что я обнаружил вторую ошибку, на этот раз в моделировании : и это связано с тем, стандартное отклонение или значение отклонения. $0.007$

В частности: Используя исправленную систему уравнений и случайное возмущение (т.е. как написано в статье ), я получаю следующий график последнего Всего 120 реализаций 3000: $\epsilon_i \sim N(0, \sigma^2 = 0.007), \implies SD = 0.0837$

Обратите внимание на значения на вертикальной оси: они намного больше, чем диапазон значений, показанный на рисунке 1 в статье авторов.

Но если я генерирую возмущения в соответствии с , то я получаю $\epsilon_i \sim N(0, \sigma^2 = 0.00049), \implies SD = 0.007$

Теперь диапазон значений совпадает с отображаемым на графике статьи. Это может быть случай, когда правильная дисперсия равна и авторы неправильно использовали ее как StDev. Но это также может быть случай, когда правильная дисперсия равна а правильное SD равно . Таким образом, моделирование было правильным (в соответствии с полученной оценкой), но по ошибке они назвали в статье «Дисперсия», что следует называть «Стандартное отклонение».$0.007$$0.000049$$0.007$

Я попытаюсь связаться с авторами по этим двум вопросам.