Linear Regression Coefficient Interpretation I have a regression problem: regress y on x1 x2
y: Sales revenue,  x1: is_holiday (dummy variable),  x2: another continuous variable (note that this is just the simplified version)
When interpreting the coefficient of x1, i.e. beta 1, we will say: "when holding all other variables constant, on a holiday, there is beta 1 dollars more revenue compared with a non-holiday.
[1] However, how does regression achieve "holding all other variables constant"? Does it achieve by taking partial derivative when solving First Order Condition of cost function (sum of squared error)?
[2] From the formula of beta 1, can we find any intuitions that the calculation is actually 'holding other variables constant'? I guess it's not observable from the formula
[3] In the dataframe, the regression is not matching holiday sample with non-holiday sample which have the same x2 value, then compare y right? I think such matching is not what regression does. I need to explain the intuition of beta 1, i.e. how "holding other variables constant" is achieved. Hopefully get some explanation in layman terms.
Thanks!
 A: The terminology "holding other variables constant" is a bit misleading, I think. I think it's more accurate to state that the coefficient $\beta_{i}$ represents the contribution of the factor $x_{i}$ that's uncorrelated to the other factors. In the following discussion, I'm going to refer to the general case of multiple quantitative variables, as I think that's more what your questions are pointing towards, instead of the two factor case with an indicator variable you use in your example. 
In order to illustrate the issue, recall that for simple linear regression the model is $y = \beta_{0} + \beta_{1}x + \varepsilon$, where $\varepsilon$ is the random error term that satisfies the usual regression assumptions. Our estimate is then $\hat{y} = \hat{\beta}_{0} + \hat{\beta}_{1}x_{1}$. In this case, as you seem to be aware, you can show that for the estimated coefficient $\hat{\beta}_{1} = \frac{\text{Cov}(y, x)}{\text{Var}(x)}$. 
Now for multiple regression, we have a model of the form $y = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + \dots+ \beta_{p}x_{p} + \varepsilon$. The estimate is of the form $\hat{y} = \hat{\beta}_{0} + \hat{\beta}_{1}x_{1} + \hat{\beta}_{2}x_{2} + \dots + \hat{\beta}_{p}x_{p}$. If you try to naively derive the same relationship for your estimated coefficients here as was found with simple linear regression, say for $x_{1}$, you find that
$\frac{\text{Cov}(y, x_{1})}{\text{Var}(x_{1})} = \frac{\text{Cov}(\hat\beta_{0} + \hat\beta_{1}x_{1}+ \hat\beta_{2}x_{2} + \dots+ \hat\beta_{p}x_{p}, x_{1})}{\text{Var}(x_{1})} 
 = \frac{\hat\beta_{1}\text{Var}(x_{1}) + \hat\beta_{2}\text{Cov}(x_{2}, x_{1})+ \dots+ \hat\beta_{p}\text{Cov}(x_{p}, x_{1})}{\text{Var}(x_{1})} \color{red}\neq \hat\beta_{1}  $
except in the special case where all the variables $x_{2}, \dots, x_{p}$ are uncorrelated to $x_{1}$. 
If we instead perform a transformation of the form (note that I have essentially set $\alpha_{1} = 0$ so that the subscript indices directly correspond, hopefully making things a little easier to understand) $x_{1} = \alpha_{0} + \alpha_{2}x_{2} + \alpha_{3}x_{3} + \dots + \alpha_{p}x_{p} + \tau$, so that $\tau$ is uncorrelated with any of $x_{2}, \dots, x_{p}$, you can perform a similar derivation to that used in simple linear regression to show that $\hat{\beta}_{1} = \frac{\text{Cov}(y, \tau)}{\text{Var}(\tau)}$. To see the relationship between the covariance of $\tau$ and $y$ with $\hat\beta_{1}$, 
$\frac{\text{Cov}(y, \tau)}{\text{Var}(\tau)} = \frac{\text{Cov}(\hat\beta_{0} + \hat\beta_{1}x_{1}+ \hat\beta_{2}x_{2} + \dots+ \hat\beta_{p}x_{p}, \tau)}{\text{Var}(\tau)} = \frac{\hat\beta_{1}\text{Cov}(x_{1}, \tau)+ \sum_{i = 2}^{p}\hat\beta_{i}\text{Cov}(x_{i}, \tau)}{\text{Var}(\tau)}$
The sum term in the numerator is equal to $0$ by our assumption that $\tau$ is uncorrelated with any of $x_{2}, \dots, x_{p}$. Substituting our transformation in for $x_{1}$ in our remaining term, we have
$\hat\beta_{1}\text{Cov}(x_{1}, \tau) = \hat\beta_{1}\text{Cov}(\alpha_{0} + \alpha_{2}x_{2} + \alpha_{3}x_{3} + \dots + \alpha_{p}x_{p} + \tau, \tau) = \hat\beta_{1}\text{Var}(\tau)$, where the last equality follows again by $\tau$ being uncorrelated with any of $x_{2}, \dots, x_{p}$. The $\text{Var}(\tau)$ term in the numerator and denominator now cancel, and we are left with $\hat\beta_{1}$, demonstrating that indeed $\hat\beta_{1} = \frac{\text{Cov}(y, \tau)}{\text{Var}(\tau)}$. But by construction, $\tau$ was the part of $x_{1}$ that was uncorrelated with any of the other variables $x_{2}, \dots, x_{p}$. The above argument also shows, mutatis mutandis, the same relationship for the other coefficients $\hat\beta_{2}, \dots, \hat\beta_{p}$. 
