Simulating Real Business Cycle
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.