Simulating Real Business Cycle Basically I need to replicate Hartley's 'A User's Guide to Solving Real Business Cycle Models' (http://www.econ.ucdavis.edu/faculty/kdsalyer/LECTURES/Ecn235a/Linearization/ugfinal.pdf). Specifically, 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.

I simulate it using the following logic: say at time $t$, everything is at steady state and all the values are 0, from which we have $k_{t+1}$. Then, at $t+1$ by giving a shock to the system through $\varepsilon$, 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:



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

n&lt;-300

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

for (ii in 1:n){

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

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

  solution &lt;- solve(lmat, rmat)
  ct1 &lt;- solution[1, ]
  ht1 &lt;- solution[2, ]

  ##now solve for yt1 and kt2 and it1
  yt1 &lt;- zt1 + 0.36 * kt1 + 0.64 * ht1
  kt2 &lt;- -0.07 * ct1 + 1.01 * kt1 + 0.06 * ht1 + 0.1 * zt1
  it1 &lt;- 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 &lt;- rbind(data.simulated, data.table(t = ii + 1, kval = kt2), fill = TRUE)
}


a &lt;- data.simulated[, list(t, cval, ival, yval)]
a &lt;- 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 specified in the paper is inherently unstable and ergo the results, or did I make a mistake somewhere?
Explosiveness

The paper contains an error, which causes the explosive dynamics in your simulation (although presumably the underlying computations in the paper were correct). The equilibrium condition derived from eigenvalue decomposition is contained in the third row of matrix $Q^{-1}$ on page 12 of the paper, with variables ordered as $(c,k,h,z)$ (I'll drop tildas, so all lowercase variables are to be understood as log-deviations). Comparing with eqn. (16) on p. 13, we see that coefficients for $k$ and $h$ are switched, and so the correct condition is

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

Simulation

First, we can express consumption and labor as linear function of state variables (no need to solve the system at each step of the simulation). The intertemporal and intratemporal equilibrium conditions can be written as

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

so after multiplying by an inverse we get

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

Next, transition for states can be written as

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

which can be reduced by substuting for control variables to

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

Now the simulation should be trivial, here's a Matlab/Octave example:

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




Of course in practice, you should probably recompute the whole solution, including the eigenvalue decomposition, so that you would be able to change parameters, etc.