Assume that the true causal relation is
$$x_i = ay_i + u_i \tag{1}$$
with the $u$-vector independent of the $y_i$-vector, but we mispecify
$$y_i = bx_i + \epsilon_i \tag{2}$$
And we get the theoretical relationship (substituting $(1)$ in $(2)$ and applying expected values)
$$b = \frac 1a \tag{3}$$
Attempting an OLS estimation for $b$ we get
$$\hat b = \frac {\sum x_iy_i}{\sum x_i^2}$$
What does this estimate in reality?
We need to plug in eq.$(1)$ to find out (since this is the true causal relationship, by assumption, while $(2)$ is just a figment of our imagination). We get
$$\hat b = \frac {\sum (ay_i + u_i)y_i}{\sum (ay_i + u_i)^2} = \frac {a\sum y_i^2 + \sum u_iy_i}{a^2\sum y_i^2 + 2a\sum u_iy_i + \sum u_i^2}$$
This is certainly a biased estimator. Asymptotically, given the independence between $y_i$ and $u_i$ (orthogonality would suffice) we will get (multiplying up and down by $(1/n)$)
$$\hat b \xrightarrow{p} \hat b_p = \frac {aE(y^2)}{a^2E(y^2) + \sigma_u^2} = \frac 1a \cdot \left(\frac {E(y^2)}{E(y^2)+ (\sigma_u/a)^2}\right) \tag{4}$$
This shows that $\hat b$ is an inconsistent estimator for $1/a$. The term in the big parenthesis is always positive and smaller than unity, so we get the "attenuation bias" (bias towards zero) phenomenon, i.e. the plim of $\hat b$ will be closer to the zero value than $1/a$ irrespective of whether $a$ is positive or negative.
Can we do anything else? Well what if we attempt to estimate the variance using the residuals? We have
$$\hat \sigma^2_{\epsilon} = \frac 1n \sum \hat \epsilon_i^2 = \frac 1n \sum [y_i-\hat b(ay_i +u_i)]^2 = \frac 1n \sum [(1-\hat b a)y_i-\hat bu_i)]^2$$
$$= (1-\hat ba)^2\frac 1n \sum y_i^2 -2(1-\hat b a)\hat b\frac 1n\sum y_iu_i + \hat b^2\frac 1n\sum u_i^2$$
The probability limit of this is
$$\hat \sigma^2_p = (1-\hat b_pa)^2E(y^2) + \hat b_p^2 \sigma^2_u \tag{5}$$
Now note that
a) for the left-hand sides of $(4)$ and $(5)$ we have consistent estimates from the estimation procedure (since they are the actual probability limits of the estimators we used)
b) We can estimate consistently $E(y^2)$
So if you rearrange $(4)$ to solve for $a$, rearrange $(5)$ to solve for $\sigma^2_u$, and use $\hat b$ instead of $\hat b_p$, $\hat \sigma^2_{\epsilon}$ instead of $\hat \sigma^2_p$, and $(1/n)\sum y_i^2$ instead of $E(y^2)$ you have a system of two equations in two unknowns ($a$ and $\sigma^2_u$).
Does it give a solution?
I am answering under the supervision of CV's peers. Be critical.
Assume one has the following model specification
$\boldsymbol{y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{u}$
where $\boldsymbol{y}$ is a $n \times 1$ vector, $\boldsymbol{X}$ an $n \times k$ matrix, $\boldsymbol{\beta}$ a $k \times 1$ hyperparameter and $\boldsymbol{u}$ a $n \times 1$ vector of homoscedastic but autocorrelated residuals. At this stage we still do not know how those are autocorrelated.
Assume that one omited to include another variable, say a $n \times 1$ vector $\boldsymbol{z}$, whose endogeneity consists of its autocorrelation such that
$\boldsymbol{z} = f(\boldsymbol{z})$
where $f$ is assumed to be a bijective/invertible vector function which specifies the correlation structure between the $n$ components $z_{i=1,...,n}$ of $\boldsymbol{z}$.
This means that $\boldsymbol{u}$ is hiddenly generated as follows (where $\gamma \neq 0$ is a scalar parameter and $\boldsymbol{v}$ is $n \times 1$ vector of errors assumed to be iid normal.)
$\boldsymbol{u} = \gamma\boldsymbol{z} + \boldsymbol{v} \iff \frac{1}{\gamma}(\boldsymbol{u}-\boldsymbol{v}) = \boldsymbol{z} \iff f(\frac{1}{\gamma}(\boldsymbol{u}-\boldsymbol{v})) = f(\boldsymbol{z})$
But since one has $\boldsymbol{z} = f(\boldsymbol{z})$, the above last equivalence can be turned into an equality. Which leads to
$\frac{\boldsymbol{u}-\boldsymbol{v}}{\gamma} = f(\frac{\boldsymbol{u}-\boldsymbol{v}}{\gamma}) \iff u = f(\frac{\boldsymbol{u}-\boldsymbol{v}}{\gamma})\gamma+\boldsymbol{v}$
Which shows that even if the correlation is not the same as the one there is between the components of $\boldsymbol{z}$, it does exist.
Thus yes serial correlation does have something to do with endogeneity when, e.g., this endogeneity consists of an omited autocorrelated variable whose autocorrelation structure is invertible.
But actually, it is very unlikely that $f$ be invertible. I mean that, if autocorrelation works through time, $f$ is the
backshift operator, and it is not invertible.
Best Answer
This is quite a well known issue in economics and I would say that you would most likely have to estimate your equation by IV rather than OLS since your error term will most likely be correlated with your regressors, i.e. $E\left(x_{i}error_{i}\right)\neq0 $. In other words you'll have an endogeneity problem due to simultaneity bias ( http://en.wikipedia.org/wiki/Endogeneity_(applied_statistics) ). Often lagged values of the regressors are useable as instruments. There are of course other options in order to estimate the relationship at hand! Note that if your coefficient is 0.60 or not does not play any role since your estimates will be biased and inconsistent in case you have a simultaneity problem which you do.
Answering your questions step-by-step:
Does this model suffer from reverse causality? Yes it does as explained above.
In other words, is it the case that the relationship is because higher consumption is driving down income? Income affects consumption and consumption affects income as is known from economic theory.
Can I use this as a rule-of-thumb to rule out the reverse causality in this case? No. This is because your estimates are inconsistent and biased.
Is this generalizable to other cases with two variables? No it is not as it really depends on each case. There is no rule of thumb in this case other than using the Hausman test to test for simultaneity bias. What we are testing is whether or not the OLS estimates are consistent or not.
Try to look at: Campbell, John Y. and Mankiw, N. Gregory. 1989. “Consumption, Income and Interest Rates: Reinterpreting the Time Series Evidence”.