Simulación del ciclo económico real

Question

Básicamente necesito replicar el libro de Hartley 'A User's Guide to Solving Real Business Cycle Models' ( http://www.econ.ucdavis.edu/faculty/kdsalyer/LECTURES/Ecn235a/Linearization/ugfinal.pdf ). En concreto, quiero simular el sistema dinámico implícito en el modelo que se especifica como sigue:

donde $c$ es el consumo, $h$ es la oferta de trabajo, $k$ es el capital, $z$ es el proceso tecnológico autorregresivo, $y$ es la salida y $i$ es la inversión.

Lo simulo utilizando la siguiente lógica: digamos que en el momento $t$ todo está en estado estacionario y todos los valores son 0, por lo que tenemos $k_{t+1}$ . Entonces, en $t+1$ dando una sacudida al sistema a través de $\varepsilon$ , resuelvo para $c_{t+1}$ y $h_{t+1}$ (ya que tengo la "conmoción $z_{t+1}$ y previamente obtenido $k_{t+1}$ . Luego, conecto esos dos para recuperar el resto, es decir $y_{t+1}, i_{t+1}, k_{t+2}$ y repetir el proceso.

Desgraciadamente, obtengo un proceso explosivo que no tiene sentido:

También incluyo el código R que se utiliza para simular esto:

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 mi pregunta es simple - es el sistema especificado en el papel es intrínsecamente inestable y ergo los resultados, o he cometido un error en alguna parte?

Answer 1

Explosividad

El documento contiene un error, que provoca la dinámica explosiva en su simulación (aunque presumiblemente los cálculos subyacentes en el documento eran correctos). La condición de equilibrio derivada de la descomposición de valores propios está contenida en la tercera fila de la matriz $Q^{-1}$ en la página 12 del documento, con las variables ordenadas como $(c,k,h,z)$ (Dejaré de lado las tildes, por lo que todas las variables en minúsculas deben entenderse como desviaciones logarítmicas). Comparando con la ecuación (16) de la página 13, vemos que los coeficientes de $k$ y $h$ se conmutan, por lo que la condición correcta es

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

Simulación

En primer lugar, podemos expresar el consumo y el trabajo como función lineal de las variables de estado (no es necesario resolver el sistema en cada paso de la simulación). Las condiciones de equilibrio intertemporal e intratemporal pueden escribirse como

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

por lo que después de multiplicar por un inverso obtenemos

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

A continuación, la transición de los estados puede escribirse como

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

que puede reducirse sustituyendo las variables de control por

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

Ahora la simulación debería ser trivial, aquí hay un ejemplo de 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')

Por supuesto, en la práctica, probablemente habría que volver a calcular toda la solución, incluida la descomposición de valores propios, para poder cambiar los parámetros, etc.

Answer 2

Noticias finales 20 de marzo de 2015 : He enviado un correo electrónico al prof. K. Salyer, uno de los autores de la Guía del Usuario. En una comunicación repetida, verificó que ambos problemas (ver mi respuesta más abajo, así como la respuesta de @ivansml), existen:

a) La ecuación correcta para la ley del movimiento del consumo es la que muestra @ivansml

b) El número $0.007$ es erróneamente llamada "varianza" (p. 11) en el documento. En realidad, es la desviación estándar y, de hecho, tal magnitud es un hallazgo típico en los datos (prof. Salyer hace referencia a la página 22 del capítulo 1 de Frontiers of Business Cycle Research de Cooley y Prescott).

Ambos errores se refieren a la versión impresa del periódico. En otras palabras, las simulaciones de la Figura 1 del documento son correctas: utilizan la ecuación correcta para el consumo, y utilizan $0.007$ como la desviación estándar de la perturbación en el proceso tecnológico. Por lo tanto, es el segundo gráfico a continuación que es válido.

---

FASE A
He comprobado por simulación (y utilizando la desviación estándar correcta) que el modelo explota, aunque lo hace hacia arriba en lugar de hacia abajo. Debe haber un error de cálculo en el documento, que sin embargo no se "transmitió" a las simulaciones de los autores. Por el momento no se me ocurre nada más, ya que la metodología es estándar. Estoy intrigado y por eso sigo trabajando en ello.

FASE B
Después de la respuesta de @ivansml, que identificó un error en el documento (que aparentemente no se cometió en las simulaciones de los autores) Creo que He identificado un segundo error, esta vez en las simulaciones y está relacionado con el hecho de que $0.007$ es una desviación estándar o un valor de varianza.

Específicamente: Utilizando el sistema de ecuaciones corregido, y una perturbación aleatoria $\epsilon_i \sim N(0, \sigma^2 = 0.007), \implies SD = 0.0837$ (es decir, como escrito en el papel ) obtengo el siguiente gráfico de las últimas 120 realizaciones de un total de 3.000 :

Observe los valores en el eje vertical: son mucho mayor que el rango de valores que aparece en la figura 1 del artículo de los autores.

Pero si genero perturbaciones de acuerdo a $\epsilon_i \sim N(0, \sigma^2 = 0.00049), \implies SD = 0.007$ , entonces obtengo

Ahora el rango de valores coincide con los que aparecen en el gráfico del documento. Puede darse el caso de que la Varianza correcta sea $0.007$ Pero también puede ser que la varianza correcta sea $0.000049$ y la SD correcta es $0.007$ . Así que la simulación era correcta (en línea con la estimación obtenida), pero por error llamaron en el documento "Varianza" lo que debería llamarse "Desviación Estándar".

Intentaré ponerme en contacto con los autores sobre estos dos asuntos.