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:

enter image description here

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:

enter image description here

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


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?


2 Answers 2



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


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()];
Y = [0.53 0.47; -0.47 1.47] * X;
plot(1:T, [X; Y])


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.

  • $\begingroup$ (+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). $\endgroup$ Commented Mar 20, 2015 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.

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.

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 : enter image description here

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 enter image description here

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.

  • $\begingroup$ 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. $\endgroup$
    – Sarunas
    Commented Mar 18, 2015 at 19:52
  • $\begingroup$ @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. $\endgroup$ Commented Mar 18, 2015 at 20:14
  • $\begingroup$ 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. $\endgroup$
    – Sarunas
    Commented Mar 18, 2015 at 20:19
  • $\begingroup$ thanks a bunch for your effort. you helped the hell out of me! at least, it was not me who made the mistake (fundamentally) :) $\endgroup$
    – Sarunas
    Commented Mar 19, 2015 at 6:29
  • 1
    $\begingroup$ @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) $\endgroup$
    – ivansml
    Commented Mar 19, 2015 at 16:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.