Simulating a dynamical system

Question

Basically I need to replicate Hartley's 'A User's Guide to Solving Real Business Cycle Models' . Specifically (to make question relevant to stats.stackexchange), I want to simulate the dynamical system implied by the model which is specified as follows:

where $c$ is consumption, $h$ is labour supply, $k$ is capital, $z$ is the autoregressive technological process, $y$ is the output and $i$ is investment. The important point is that these represent percentage deviations from steady state (and so are growth rates) and shocks are initiated through $z_t$ - autoregressive process.

I simulate it using the following logic: say at time $t$, everything is at steady state and all the values are 0 (there are no deviations from steady state), from which we have $k_{t+1}$. Then, at $t+1$ by giving a shock to the system through $\varepsilon$ (which is assumed to be normally distributed with mean 0 and sd 0.007), i solve for $c_{t+1}$ and $h_{t+1}$ (as I have the 'shocked' $z_{t+1}$ and previously obtained $k_{t+1}$. Then, I plug those two to retrieve the rest, namely - $y_{t+1}, i_{t+1}, k_{t+2}$ and repeat the process.

Unfortunately, I get an explosive process which doesn't make sense (the series should be stationary as implied by the economic theory):

I also include R code that is used to simulate this:

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 my question is simple - is the system that is specified by the paper is inherently unstable? I'm not sure how could check it analytically, so hopefully someone could help.

Answer 1

The way you have proceeded sounds reasonable to me. Another approach to do it is like this.

You have your system of equations (I have rearranged the order the variables appear in each equation and removed the tildes) $$ c_t=0.02k_t+0.54h_t+0.44z_t\\ k_{t+1}=-0.07c_t+1.01k_t+0.06h_t+0.10z_t\\ h_t=-2.78c_t+k_t+2.78z_t\\ y_t=0.36k_t+0.64h_t+z_t\\ i_t=-2.92c_t+3.92y_t\\ z_t=0.95z_{t-1}+\epsilon_t. $$ Write this using matrices: $$ \begin{pmatrix}1& -0.02 & -0.54 & 0 & 0 &-0.44\\ 0 & 1 & 0 & 0 & 0& 0\\ 2.78&-1&1&0&0&-2.78\\ 0&-0.36&-0.64&1&0&-1\\ 2.92&0&0&-3.92&1&0\\ 0& 0&0&0&0&1\end{pmatrix}\begin{pmatrix}c_{t+1} \\k_{t+1}\\h_{t+1}\\y_{t+1}\\i_{t+1}\\z_{t+1}\\ \end{pmatrix}=\begin{pmatrix}0& 0 & 0 & 0 & 0 &0\\ -0.07 & 1.01 & 0.06 & 0 & 0& 0.10\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0.95\end{pmatrix}\begin{pmatrix}c_{t} \\k_{t}\\h_{t}\\y_{t}\\i_{t}\\z_{t}\end{pmatrix}+\begin{pmatrix}0\\0\\0\\0\\0\\\epsilon_{t+1}\end{pmatrix} $$

In vector notation, your model is $$ Ax_{t+1}=Bx_t+e_t $$ so premultiplying by the inverse of $A$ $$ x_{t+1}=A^{-1}Bx_t+A^{-1}e_t $$

We might as well simulate the system from this relation. However, here we can note that the largest eigenvalue of $A^{-1}B$ is 1.01. This means that the system is explosive. So the (allegedly) strange behavior you have noticed does not seem to be a result of your implementation, but of the system itself.

A <- matrix(c(1, -0.02, -0.54, 0, 0, -0.44,
              0, 1, 0, 0, 0, 0,
              2.78, -1, 1, 0, 0, -2.78,
              0, -0.36, -0.64, 1, 0, -1,
              2.92, 0, 0, -3.92, 1, 0,
              0, 0, 0, 0, 0, 1), ncol = 6, byrow = TRUE)

B <- matrix(c(rep(0, 6),
              -0.07, 1.01, 0.06, 0, 0, 0.10,
              rep(0, 6*3),
              0, 0, 0, 0, 0, 0.95), ncol = 6, byrow = TRUE)

## Largest eigenvalue is > 1, so system is explosive
eigen(solve(A) %*% B)$values

n <- 300
x <- matrix(0, nrow = 6, ncol = n)
for (ii in 2:n) {
    e <- c(rep(0, 5), rnorm(1, 0, 0.007))
    x[, ii] <- solve(A) %*% B %*% x[, ii-1] + solve(A) %*% e
}