Symulowanie prawdziwego cyklu koniunkturalnego

Zasadniczo muszę powielić „Przewodnik użytkownika dotyczący rozwiązywania rzeczywistych modeli cyklu koniunkturalnego” ( http://www.econ.ucdavis.edu/faculty/kdsalyer/LECTURES/Ecn235a/Linearization/ugfinal.pdf ). W szczególności chcę symulować system dynamiczny implikowany przez model, który jest określony w następujący sposób:

gdzie to konsumpcja, to podaż pracy, to kapitał, to autoregresywny proces technologiczny, to produkcja, a to inwestycja.h k z y $c$$h$$k$$z$$y$$i$

Symuluję to za pomocą następującej logiki: powiedzmy, że w czasie wszystko jest w stanie ustalonym, a wszystkie wartości to 0, z których mamy . Następnie, w czasie , powodując szok dla systemu poprzez , rozwiązuję dla i (ponieważ mam „zszokowany” i wcześniej otrzymałem . Następnie podłączam te dwa, aby odzyskać resztę, a mianowicie - i powtórzę proces.k t + 1 t + 1 ε c t + 1 h t + 1 z t + 1 k t + 1 y t + 1 , i t + 1 , k t + $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}$

Niestety dostaję wybuchowy proces, który nie ma sensu:

Podaję również kod R, który służy do symulacji tego:

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()

Sy moje pytanie jest proste - czy system określony w artykule jest z natury niestabilny i powoduje ergo wyniki, czy też gdzieś popełniłem błąd?

Wybuchowość

Artykuł zawiera błąd, który powoduje wybuchową dynamikę w twojej symulacji (chociaż przypuszczalnie podstawowe obliczenia w pracy były prawidłowe). Warunek równowagi uzyskany z rozkładu wartości własnej jest zawarty w trzecim rzędzie macierzy na stronie 12 artykułu, ze zmiennymi uporządkowanymi jako (zrzucę tyldy, więc wszystkie zmienne pisane małymi literami należy rozumieć jako odchylenia logarytmiczne). W porównaniu z eqn. (16) na str. 13 widzimy, że współczynniki i są przełączane, a więc prawidłowy warunek to$Q^{-1}$$(c,k,h,z)$$k$$h$

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

Symulacja

Po pierwsze, możemy wyrazić zużycie i pracę jako liniową funkcję zmiennych stanu (nie trzeba rozwiązywać systemu na każdym etapie symulacji). Międzyokresowe i międzyporowe warunki równowagi można zapisać jako

$$\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}$$

więc po pomnożeniu przez odwrotność otrzymujemy

$$\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}$$

Następnie przejście dla stanów można zapisać jako

$$\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}$$

które można zmniejszyć, zastępując zmienne kontrolne wartością

$$\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}$$

Teraz symulacja powinna być trywialna, oto przykład 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')

Oczywiście w praktyce powinieneś prawdopodobnie ponownie obliczyć całe rozwiązanie, w tym rozkład wartości własnej, abyś mógł zmienić parametry itp.

Ostatnie wiadomości 20 marca 2015 : Wysłałem e-mail do prof. K. Salyer, jeden z autorów Przewodnika użytkownika. W powtarzającym się komunikacie potwierdził, że oba problemy (zobacz moją odpowiedź poniżej oraz odpowiedź @ivansml) istnieją:

a) Prawidłowe równanie dla prawa ruchu konsumpcji jest takie, jak pokazuje @ivansml

b) Liczba jest błędnie nazywana „wariancją” (s. 11) w pracy. W rzeczywistości jest to odchylenie standardowe i rzeczywiście taka wielkość jest typowym odkryciem w danych (prof. Salyer przywołał „s. 22 z rozdz. 1 Cooley i Prescott's Frontiers of Business Cycle Research).$0.007$

Oba błędy dotyczą wydrukowanej wersji artykułu. Innymi słowy, symulacje przedstawione na rycinie 1 artykułu są poprawne: wykorzystują prawidłowe równanie zużycia i wykorzystują jako odchylenie standardowe zakłócenia w procesie technologicznym. Więc jest to drugi wykres poniżej, że jest poprawny.$0.007$

---

FAZA A Za
pomocą symulacji (i przy użyciu prawidłowego odchylenia standardowego) zweryfikowałem, że model eksploduje, chociaż robi to raczej w górę niż w dół. W artykule musi być błąd obliczeniowy, który jednak nie został „przesłany” do symulacji autorów. W tej chwili nie mogę wymyślić nic innego, ponieważ metodologia jest standardowa. Jestem zaintrygowany i nadal nad tym pracuję.

FAZA B
Po odpowiedzi @ ivansml, która zidentyfikowała błąd w pracy (który najwyraźniej nie został popełniony w symulacjach autorów) , myślę , że zidentyfikowałem drugi błąd, tym razem w symulacjach : i jest to związane z tym, czy jest odchyleniem standardowym lub wartością wariancji. $0.007$

W szczególności: używając skorygowanego układu równań i przypadkowego zaburzenia (tj. Jak napisano w pracy ) otrzymuję następujący wykres ostatniego 120 realizacji na 3000 ogółem: $\epsilon_i \sim N(0, \sigma^2 = 0.007), \implies SD = 0.0837$

Zwróć uwagę na wartości na osi pionowej: są one znacznie większe niż zakres wartości pokazany na rycinie 1 w pracy autora.

Ale jeśli zakłócenia zgodnie z , to otrzymuję $\epsilon_i \sim N(0, \sigma^2 = 0.00049), \implies SD = 0.007$

Teraz zakres wartości odpowiada wartościom pojawiającym się na wykresie papieru. Może się zdarzyć, że poprawna wariancja wynosi a autorzy nieprawidłowo zastosowali ją jako StDev, ale może się też zdarzyć, że poprawna wariancja to a poprawna SD to . Tak więc symulacja była poprawna (zgodnie z uzyskanym oszacowaniem), ale przez pomyłkę nazwali w pracy „Wariancją”, co należy nazwać „odchyleniem standardowym”.0,000049 $0.007$$0.000049$$0.007$

Spróbuję skontaktować się z autorami w tych dwóch sprawach.