The first sentence of this wiki page claims that "In econometrics, an endogeneity problem occurs when an explanatory variable is correlated with the error term.1 "
My question is that how can this ever happen? Isn't regression beta chosen such that the error term is orthogonal to the column space of the design matrix?
You are conflating two types of "error" term. Wikipedia actually has an article devoted to this distinction between errors and residuals.
In an OLS regression, the residuals (your estimates of the error or disturbance term) $\hat \varepsilon$ are indeed guaranteed to be uncorrelated with the predictor variables, assuming the regression contains an intercept term.
But the "true" errors $\varepsilon$ may well be correlated with them, and this is what counts as endogeneity.
To keep things simple, consider the regression model (you might see this described as the underlying "data generating process" or "DGP", the theoretical model that we assume to generate the value of $y$):
$$y_i = \beta_1 + \beta_2 x_i + \varepsilon_i$$
There is no reason, in principle, why $x$ can't be correlated with $\varepsilon$ in our model, however much we would prefer it not to breach the standard OLS assumptions in this way. For example, it might be that $y$ depends on another variable that has been omitted from our model, and this has been incorporated into the disturbance term (the $\varepsilon$ is where we lump in all the things other than $x$ that affect $y$). If this omitted variable is also correlated with $x$, then $\varepsilon$ will in turn be correlated with $x$ and we have endogeneity (in particular, omitted-variable bias).
When you estimate your regression model on the available data, we get
$$y_i = \hat \beta_1 + \hat \beta_2 x_i + \hat \varepsilon_i$$
Because of the way OLS works*, the residuals $\hat \varepsilon$ will be uncorrelated with $x$. But that doesn't mean we have avoided endogeneity — it just means that we can't detect it by analysing the correlation between $\hat \varepsilon$ and $x$, which will be (up to numerical error) zero. And because the OLS assumptions have been breached, we are no longer guaranteed the nice properties, such as unbiasedness, we enjoy so much about OLS. Our estimate $\hat \beta_2$ will be biased.
$(*)$ The fact that $\hat \varepsilon$ is uncorrelated with $x$ follows immediately from the "normal equations" we use to choose our best estimates for the coefficients.
If you are not used to the matrix setting, and I stick to the bivariate model used in my example above, then the sum of squared residuals is $S(b_1, b_2) = \sum_{i=1}^n \varepsilon_i^2 = \sum_{i=1}^n (y_i-b_1 - b_2 x_i)^2$ and to find the optimal $b_1 = \hat \beta_1$ and $b_2 = \hat \beta_2$ that minimise this we find the normal equations, firstly the first-order condition for the estimated intercept:
$$\frac{\partial S}{\partial b_1} = \sum_{i=1}^n -2(y_i-b_1 - b_2 x_i) = -2 \sum_{i=1}^n \hat \varepsilon_i = 0$$
which shows that the sum (and hence mean) of the residuals is zero, so the formula for the covariance between $\hat \varepsilon$ and any variable $x$ then reduces to $\frac{1}{n-1} \sum_{i=1}^n x_i \hat \varepsilon_i$. We see this is zero by considering the first-order condition for the estimated slope, which is that
$$\frac{\partial S}{\partial b_2} = \sum_{i=1}^n -2 x_i (y_i-b_1 - b_2 x_i) = -2 \sum_{i=1}^n x_i \hat \varepsilon_i = 0$$
If you are used to working with matrices, we can generalise this to multiple regression by defining $S(b) = \varepsilon' \varepsilon = (y-Xb)'(y-Xb)$; the first-order condition to minimise $S(b)$ at optimal $b = \hat \beta$ is:
$$\frac{dS}{db}(\hat\beta) = \frac{d}{db}\bigg(y'y - b'X'y - y'Xb + b'X'Xb\bigg)\bigg|_{b=\hat\beta} = -2X'y + 2X'X\hat\beta = -2X'(y - X\hat\beta) = -2X'\hat \varepsilon = 0$$
This implies each row of $X'$, and hence each column of $X$, is orthogonal to $\hat \varepsilon$. Then if the design matrix $X$ has a column of ones (which happens if your model has an intercept term), we must have $\sum_{i=1}^n \hat \varepsilon_i = 0$ so the residuals have zero sum and zero mean. The covariance between $\hat \varepsilon$ and any variable $x$ is again $\frac{1}{n-1} \sum_{i=1}^n x_i \hat \varepsilon_i$ and for any variable $x$ included in our model we know this sum is zero, because $\hat \varepsilon$ is orthogonal to every column of the design matrix. Hence there is zero covariance, and zero correlation, between $\hat \varepsilon$ and any predictor variable $x$.
If you prefer a more geometric view of things, our desire that $\hat y$ lies as close as possible to $y$ in a Pythagorean kind of way, and the fact that $\hat y$ is constrained to the column space of the design matrix $X$, dictate that $\hat y$ should be the orthogonal projection of the observed $y$ onto that column space. Hence the vector of residuals $\hat \varepsilon = y - \hat y$ is orthogonal to every column of $X$, including the vector of ones $\mathbf{1_n}$ if an intercept term is included in the model. As before, this implies the sum of residuals is zero, whence the residual vector's orthogonality with the other columns of $X$ ensures it is uncorrelated with each of those predictors.
But nothing we have done here says anything about the true errors $\varepsilon$. Assuming there is an intercept term in our model, the residuals $\hat \varepsilon$ are only uncorrelated with $x$ as a mathematical consequence of the manner in which we chose to estimate regression coefficients $\hat \beta$. The way we selected our $\hat \beta$ affects our predicted values $\hat y$ and hence our residuals $\hat \varepsilon = y - \hat y$. If we choose $\hat \beta$ by OLS, we must solve the normal equations and these enforce that our estimated residuals $\hat \varepsilon$ are uncorrelated with $x$. Our choice of $\hat \beta$ affects $\hat y$ but not $\mathbb{E}(y)$ and hence imposes no conditions on the true errors $\varepsilon = y - \mathbb{E}(y)$. It would be a mistake to think that $\hat \varepsilon$ has somehow "inherited" its uncorrelatedness with $x$ from the OLS assumption that $\varepsilon$ should be uncorrelated with $x$. The uncorrelatedness arises from the normal equations.
The data generating process is:
$$ y_i = a + b_1x_{i,1} + b_2x_{i,2} + \epsilon_i$$
If we ran that regression, we would get estimates $\hat{a}$, $\hat{b}_1$, and $\hat{b}_2$, and with enough data, they would converge on $a$, $b_1$, and $b_2$ respectively.
(Technical note: We need a little randomness so we don't buy exactly one bun for each burger we buy at every visit to the grocery store. If we did this, $x_1$ and $x_2$ would be collinear.)
Now let's consider the model:
$$ y_i = a + b_1x_{i,1} + u_i $$
Observe that $u_i = b_2x_{i,2} + \epsilon_i$. Hence \begin{align*} \operatorname{Cov}(x_{1}, u) &= \operatorname{Cov}(x_1,b_2x_2 + \epsilon )\\ &= b_2 \operatorname{Cov}(x_{1},x_2) + \operatorname{Cov}(x_{1},\epsilon) \\ &= b_2 \operatorname{Cov}(x_{1},x_2) \end{align*}
Is this zero? Almost certainly not! The purchase of burgers $x_1$ and the purchase of buns $x_2$ are almost certainly correlated! Hence $u$ and $x_1$ are correlated!
If you tried to run:
$$ y_i = \hat{a} + \hat{b}_1 x_{i,1} + \hat{u}_i $$
Your estimate $\hat{b}_1$ would almost certainly be a poor estimate of $b_1$ because the OLS regression estimates $\hat{a}, \hat{b}, \hat{u}$ would be constructed so that $\hat{u}$ and $x_1$ are uncorrelated in your sample. But the actual $u$ is correlated with $x_1$ in the population!
What would happen in practice if you did this? Your estimate $\hat{b}_1$ of the price of burgers would ALSO pickup the price of buns. Let's say every time you bought a \$1 burger you tended to buy a \$0.50 bun (but not all the time). Your estimate of the price of burgers might be \$1.40. You'd be picking up the burger channel and the bun channel in your estimate of the burger price.
Suppose that we're building a regression of the weight of an animal on its height. Clearly, the weight of a dolphin would be measured differently (in different procedure and using different instruments) from the weight of an elephant or a snake. This means that the model errors will be dependent on the height, i.e. explanatory variable. They could be dependent in many different ways. For instance, maybe we tend to slightly overestimate the elephant weights and slightly underestimate the snake's, etc.
So, here we established that it is easy to end up with a situation when the errors are correlated with the explanatory variables. Now, if we ignore this and proceed to regression as usual, we'll notice that the regression residuals are not correlated with the design matrix. This is because, by design the regression forces the residuals to be uncorrelated. Note, also that residuals are not the errors, they're the estimates of errors. So, regardless of whether the errors themselves are correlated or not with the independent variables the error estimates (residuals) will be uncorrelated by the construction of the regression equation solution.
