Since this is clearly a homework problem, I won't give away the full answer. First, I'll note that the formula you posted does not appear to be correct; the denominator should be $\sum(D_i-\bar D)^2$. $\sum(D_i-\bar D)$ is just $0$.
The formula says that a sample estimate of $\beta_1$, namely, $\hat \beta_1$, is equal to the true value of $\beta_1$ in the population plus a bias term. The bias term is a ratio of two quantities, which involve $D_i$, the binary predictor; its mean $\bar D$ (also equal to the proportion of $1$s; and $\epsilon_i$, a random disturbance. The bias term can be written in terms of population moments (e.g., mean, variance, covariance, skew, etc.), and that is what you are being asked to figure out how to do.
I'll provide two insights that might help you:
- $\epsilon_i$ is centered at $0$, meaning $\epsilon_i = \epsilon_i - \bar \epsilon$.
- You can multiply the numerator and denominator by $\frac{1}{n}$, that is, $\frac{1/n}{1/n} = 1$, which seems trivial but is actually critical in converting the numerator and denominator into values that resemble population parameters.
Try using these clues to solve the problem. Let me know if you need more hints.
The original formula you presented was
$$\hat\beta_1=\beta_1 + \frac{\sum(D_i - \bar D)\epsilon_i}{\sum (D_i - \bar D)^2}$$
When you replace $\epsilon_i$ with $\epsilon_i - \bar \epsilon$, the numerator is $\sum(D_i - \bar D)(\epsilon_i - \bar \epsilon)$. Hopefully that looks like a familiar formula, though, again, you need the second insight to make the connection. Using the second insight, you can write
$\frac{\sum W}{\sum Z}$ as $\frac{\frac{1}{n}\sum W}{\frac{1}{n}\sum Z}$ (here I simply used $W$ and $Z$ as placeholders). So you can turn a ratio of sums into a ratio of means. What are they the means of in this case?
Here's what you are supposed to see:
$$
\sum(D_i - \bar D)\epsilon_i = \sum(D_i - \bar D)(\epsilon_i - \bar \epsilon)
$$
$$\begin{aligned}
\hat\beta_1 &= \beta_1 + \frac{\sum(D_i - \bar D)\epsilon_i}{\sum (D_i - \bar D)^2} \\
&= \beta_1 + \frac{\sum(D_i - \bar D)(\epsilon_i - \bar \epsilon)}{\sum (D_i - \bar D)^2} \\
&= \beta_1 + \frac{\frac{1}{n}\sum(D_i - \bar D)(\epsilon_i - \bar \epsilon)}{\frac{1}{n}\sum (D_i - \bar D)^2} \\
&= \beta_1 + \frac{\text{Cov}(D, \epsilon)}{\text{Var}(D)}
\end{aligned}$$
That's it. $\text{Cov}(D, \epsilon)$ and $\text{Var}(D)$ are moments of the joint distribution of $D$ and $\epsilon$. The usual OLS estimator for univariable regression is $\frac{\text{Cov}(D, Y)}{\text{Var}(D)}$; expanding $Y$ into $Y = \beta_0 + \beta_1 D + \epsilon$ will reduce to the formula I derived above.
Expanding the equation for $Y$ into the usual OLS solution:
$$\begin{aligned}
\hat\beta_1 &= \frac{\text{Cov}(D, Y)}{\text{Var}(D)} \\
& = \frac{\text{Cov}(D, \beta_0 + \beta_1 D + \epsilon)}{\text{Var}(D)} \\
& = \frac{\text{Cov}(D, \beta_0) + \text{Cov}(D, \beta_1 D) + \text{Cov}(D, \epsilon)}{\text{Var}(D)} \\
&= \frac{0 + \beta_1 \text{Cov}(D, D) + \text{Cov}(D, \epsilon)}{\text{Var}(D)} \\
&= \frac{\beta_1 \text{Var}(D)}{\text{Var}(D)} + \frac{\text{Cov}(D, \epsilon)}{\text{Var}(D)} \\
&= \beta_1 + \frac{\text{Cov}(D, \epsilon)}{\text{Var}(D)}
\end{aligned}$$