Why does the parameter variance change when control variables are added to a regression model? If I add a control variable to my regression, this changes the variances of the parameter estimates. Why is this the case? Is this because SSE(=explained sum of squares) increases and therefore the SSR(=residual sum of squares) decreases?  So basically as long as I avoid multicollinearity adding more control variables will always improve the statistical significance of the parameter estimates?
 A: The question is what happens to the variance of a parameter estimates ($\hat \beta_i$) with the introduction of a new control variable in a regression model? 
The variance of the estimated parameter will increase.
The more regressors we add, the lower the $\text{RSS}$ (residual sum of squares) and the higher the $R^2$. In fact, if we just add columns of pure noise (rnorm()) to our data, and run OLS regressions after every addition, we will see a marked, monotonic decrease in the $\text{RSS}$. I tested this with an example:

The idea is that in OLS the hat matrix $X (X^TX)^{-1}X^T$ is a projection matrix of the vector of observed $y$ values onto the column space of the model matrix. The higher the dimensions of the column space (the more vectors to form a basis), the closer $\hat y$ will be to $y$. But at a price!
The estimation of the variance of $\hat \beta_i$ is given by:
$Var[\hat\beta_i]= \sigma^2(X^TX)^{-1}$ with $\sigma^2$ corresponding to the variation of the observations around the predicted values, ($Var[\epsilon|X] = \sigma^2I_n$). The estimation of $\sigma^2$ from the sample is $s^2=\frac{e^Te}{n-p}$ or $\text{MSE}$ (mean squared error). The denominator is the number of observations, $n$ minus the number of parameters, $p$, counting the intercept. It is also referred to as the error or residual degrees of freedom ($\small\text{no. observations−no. indepen't variables−1}$). Alternatively, the formula for the $\sigma^2$ estimation, can be expressed as $\text{MSE = RSS/df}$ with $\text{RSS}$ being the same as $\text{SSE}$ (sum of squared errors), $\sum_1^n(y_i - \hat y)^2$. 
Therefore the estimation of the variance of the parameter $\hat\beta_i$ is $s^2(X^TX)^{-1}$ or $\text{MSE}\times(X^TX)^{-1}$.
And I think this may be a source of confusion - just because adding another regressor decreases the $\text{MSE}$ by decreasing the $\text{RSS}$ - if it actually does at all, because of the change in the degrees of freedom in the denominator, it can't be said that the variance for the estimated parameter decreases.
In a parallel OLS simulation (here), I found, quite anecdotally, a $\small 3.5\%$ increase in $\text{MSE}$ when adding an extra regressor, with a minimal decrease in $\text{RSS}$. 
What controlled the increased variance for the estimate of $\hat\beta_i$ turned out to be the entry for $\hat\beta_i$ in the matrix of cofactors involved in the calculation of the inverse of $(X^TX)^{-1}$, explaining a variance for the estimate of the parameter $\hat\beta_i$ $2.8$ times higher in the presence of the control variable.
An alternative formula for the variance not applicable to the intercept is:
$\large \text{var}[\hat\beta_i]= \frac{\sigma^2}{n \times var[X_i]\times(1\,-\,R_i^2)}$. The key here is that $R_i^2$ is the $R$-square of the regression of the corresponding variable for the parameter $\hat\beta_i$, or $X_i$ against all the other regressors. Therefore, the better the regression model of $X_i \sim \text{control variable}$, the higher the estimated $\text{var}[\hat\beta_i].$
A: The variances of the parameter estimates change when one adds more explanatory variables to a regression model because this action affects the overall model fit. To this why this is the case remember that:


*

*You are concerned about the maximum likelihood estimates (MLE) of your linear model.

*The Hessian of the log-likelihood evaluated at the ML estimates is the observed Fisher information.

*The estimated standard errors are the square roots of the diagonal elements of the inverse of the observed Fisher information matrix.

*ML estimates and least squares (LS) estimates are identical given you have a response variable coming from an exponential family.


OK, let's say you have already $m_1$ model such that $y = \beta_0  + \beta_{1} x_{1} + \epsilon$ and you  augment it with an additional explanatory variable $x_{2}$. Let's also assume that $x_{2}$ is not orthogonal to $y$ or $x_{1}$, (ie. $y \not\perp x_{2}$ and $x_{1} \not\perp x_{2}$ ). This is crucial because it means that when adding $x_{2}$ in the model $m_1$ we added some additional information about $y$. This increased amount of information will cause the sum-of-squares-error (SSE) to get smaller. (Geometrically speaking by incorporating $x_2$ in $m_1$ we add an additional non-orthogonal plane to project the sample $y$, thus the projections $\hat{y}$ are closer to  original data $y$ than before; see this excellent answer by silverfish here for more details.)
Let's interpreter getting a smaller SSE means in terms of LS and then ML estimates: A smaller SSE means that the optimum (the vector of $\hat{\beta}$'s) of our LS function is different than before. This means that ML estimates are different (point 4). Which means that unless the Hessian of the log-likelihood function at this new solution (point 2) are exactly the same as before, the standard errors will be different too (point 3). Therefore adding a new variable $x_2$ is extremely unlikely to result the same standard errors.
Your second assessment about improved statistical significance is directly related to the above. The wording: "As long as I avoid multicollinearity adding more control variables will always improve the statistical significance of the parameter estimates." is slightly over-optimistic. The significance of a predictor will always get smaller if you have another predictor that is collinear with it. You can also say that the significance of the others predictors will never get smaller if the newly predictor is not collinear.
Nevertheless, this last statement downplays how hard is to avoid even the slightest collinearity between two predictors while at the same time making sure that they are related to $y$. That could possibly be the case if one used PCA scores as predictors variables but this really escapes the scope of this answer. 
OK, how about a code example in R. We will define a model m1 and then we will add we will add a variable x2 to it and see if this all holds.
set.seed(1234)
N = 1000;
x1 = runif(N);
e = rnorm(N);
y = 3*x1 + e
m1 <- lm(y~x1) 

sum(m1$residuals^2)
#  898.0496

x2 = runif(N) 
cor(x2,y)  # 0.02983 linear correlation 
cor(x2,x1) # 0.02598 linear correlation 

m2 <- lm(y~x1+x2) 
sum(m2$residuals^2) 
#   897.7952

So one immediately sees that even adding a variable $x_2$ that is pure noise, exactly because of some small spurious correlations he got a smaller SSE (which was what @whuber's comment essentially was about). Clearly the improvement in the SSE is marginally (after all we added just noise) but it is still there. What about the standard errors? Are they the same?
summary(lm1)$coefficients["x1", 'Std. Error']    # [1] 0.103062  
summary(lm2)$coefficients["x1", 'Std. Error']    # [1] 0.103134
summary(lm1)$coefficients["(Intercept)", 'Std. Error']    # [1] 0.0602753
summary(lm2)$coefficients["(Intercept)", 'Std. Error']    # [1] 0.0783329

Again, while the difference is marginal, the difference is there. The standard errors for $\hat{\beta_0}$ ((Intercept)) and  $\hat{\beta_1}$  (x1) are different between the two models and that is because as the location of the optimum changed almost certainly the Hessian of the likelihood at the new optimum point is also different. In this example we also see that the standard errors are larger. This was not certain but it was likely. What about the statistical significance though?
summary(lm2)$coefficients["x1", 'Pr(>|t|)'] # [1] 9.3e-128
summary(lm1)$coefficients["x1", 'Pr(>|t|)'] # [1] 5.1e-128

summary(lm1)$coefficients["(Intercept)", 'Pr(>|t|)'] # [1] 0.4606
summary(lm2)$coefficients["(Intercept)", 'Pr(>|t|)'] # [1] 0.8192

Well the statistical significance in terms of $\beta_1$ is the same (differences of $3e-128$ are inconceivable). The statistical significance of the intercept $\beta_0$ is less and this exactly because the expected value of $x_2$ ($E\{x_2\}$) is not exactly $0$. As such some tiny variation around the intercept was explained because of $x_2$. If we had centred $x_2$ we would have the same or even higher significance for $\beta_0$ because the actual estimate of $\beta_0$ would be somewhat higher (assuming it's SE did not change proportionally too.)
