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

so after multiplying by an inverse we get

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

Next, transition for states can be written as

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

which can be reduced by substuting for control variables to

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

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.

• (+1). Perhaps it would be helpful to graph output and investment, which are usually also in the focus of interest (and contribute in validating the model, when the investment series exhibits larger variability than the output series). – Alecos Papadopoulos Mar 20 '15 at 1:44

Final NEWS March 20, 2015 : I have e-mailed prof. K. Salyer, one of the authors of the User Guide. In a repeated communication, he verified that both issues (see my answer below, as well as @ivansml answer), do exist:

a) The correct equation for the law of motion of consumption is as @ivansml shows

b) The number $0.007$ is wrongly called "variance" (p. 11) in the paper. In reality, it is the standard deviation, and indeed such a magnitude is a typical finding in the data (prof. Salyer referenced "p. 22 form ch.1 of Cooley and Prescott's Frontiers of Business Cycle Research).

Both mistakes relate to the printed version of the paper. In other words, the simulations behind Figure 1 of the paper are correct: they use the correct equation for consumption, and they use $0.007$ as the standard deviation of the disturbance in the technology process. So it is the second graph below that it is valid.

PHASE A
I verified by simulation (and using the correct standard deviation) that the model explodes, although it does so upwards rather than downwards. There must be a calculating mistake in the paper, which nevertheless was somehow not "transmitted" to the authors' simulations. For the moment I cannot think of anything else, since the methodology is standard. I am intrigued and so still working on it.

PHASE B
After @ivansml's answer, which identified a mistake in the paper (that apparently was not made in the simulations of the authors), I think I have identified a second mistake, this time in the simulations: and it is related to whether $0.007$ is a standard deviation or a variance value.

Specifically: Using the corrected system of equations, and a random disturbance $\epsilon_i \sim N(0, \sigma^2 = 0.007), \implies SD = 0.0837$ (i.e. as written in the paper) I obtain the following graph of the last 120 realizations of 3,000 total :

Note the values on the vertical axis: they are much greater than the range of values that appears in Figure 1 in the authors' paper.

But if I generate disturbances according to $\epsilon_i \sim N(0, \sigma^2 = 0.00049), \implies SD = 0.007$, then I obtain

Now the range of values matches those appearing in the paper's graph. It may be the case that the correct Variance is $0.007$ and the authors incorrectly used it as StDev.But it also may be the case that the correct variance is $0.000049$ and the correct SD is $0.007$. So the simulation was correct (in line with the obtained estimate), but by mistake they called in the paper "Variance" what should be called "Standard Deviation".

I will attempt to contact the authors on these two matters.

• i managed to figure out that the process is in fact explosive as you pointed out. i did make mistake about the variance, bus since sd is 0.083 that means even bigger variation than I initially used and the process explodes much faster. what I don't get how the author managed to simulate (as he writes) 3000 observations and provide the plot of stationary series (at the end of the paper) while the process does not exhibit this property. – Sarunas Mar 18 '15 at 19:52
• @Sarunas Check your code as follows: calculate manually the first two-three values of the various processes, using the actually generated shocks, and compare with the corresponding values that the code gives you. – Alecos Papadopoulos Mar 18 '15 at 20:14
• i did that. tried going on manually. what would be useful to know from more experienced researchers is that why would the capital process be explosive? wouldn't we want it to be stationary? how otherwise steady state could be achieved? i checked eigenvalues of the system and as you pointed out earlier - the system is in fact explosive so there is nothing wrong in the code itself. either the mistakes are in the paper or I don't understand the logic. – Sarunas Mar 18 '15 at 20:19
• thanks a bunch for your effort. you helped the hell out of me! at least, it was not me who made the mistake (fundamentally) :) – Sarunas Mar 19 '15 at 6:29
• @AlecosPapadopoulos I think coefficient on capital in linearized resource constraint can exceed one (in fact it should be equal to $1/\beta$ in this model) - what matters is the stability of $(k_t, z_t)$ process once all equilibrium relationships have been substituted in. If I take the matrices $A,B$ from the paper and plug them into Paul Klein's solab solver, I get a stable solution, so I'd say there's just some numerical typo in the paper. (if you do this yourself, beware of different notation and variable order in Klein's code) – ivansml Mar 19 '15 at 16:38