I know that the variance of a linear combination of correlated random variables can be generalized (as in Variance of linear combinations of correlated random variables). My question has to do with the variance of an additional linear combination of two or more such linear combinations. This is a different question from the one linked above because applying the variance-covariance matrix for a single linear combination is trivial, but in this extended case, you would need to apply the same variance-covariance matrix twice to difference strings of constants, so as a practical matter, being able to apply the terms in the correct order is essential.
Say you have a regression model with an output variance-covariance matrix. The variance for a given prediction would be $$ \operatorname{var}(Y_1) = \sum_{i=1}^n a_i^2\operatorname{var}(X_i) + 2\sum_{i=1}^n \sum_{j\colon j > i}^n a_ia_j\operatorname{cov}(X_i,X_j) $$ for all $a_i$ assigned to each random variable $X_i$. So for the regression, $a$ would be a vector of specific values to multiply by the vector of model coefficients $X$ (the $\beta$s) to sum up and get your prediction $Y_1$.
Now suppose you want to compare this prediction to a second prediction, $Y_2$, that comes from the same regression model, but uses a different collection of values, $b_i$. For example, if you wanted to compare individuals with a certain profile of covariate attributes against individuals with a different profile of covariate attributes. The variance for this new prediction would be $$ \operatorname{var}(Y_2) = \sum_{i=1}^n b_i^2\operatorname{var}(X_i) + 2\sum_{i=1}^n \sum_{j\colon j > i}^n b_ib_j\operatorname{cov}(X_i,X_j). $$
It is clear that the prediction of the linear combination $Y_1 - Y_2$ is simply an extended linear combination of the values for the first profile, the negative values for the second profile, and the model coefficients (repeated for the second profile). But as I'm doing the math, it appears that the variance of $Y_1 - Y_2$ is also simply a linear combination of $\operatorname{var}(Y_1)$ and the negative of $\operatorname{var}(Y_2)$ $$ \operatorname{var}(Y_1 - Y_2) = \sum_{i=1}^n a_i^2\operatorname{var}(X_i) + 2\sum_{i=1}^n \sum_{j\colon j > i}^n a_ia_j\operatorname{cov}(X_i,X_j) + \sum_{i=1}^n (-b_i)^2\operatorname{var}(X_i) + 2\sum_{i=1}^n \sum_{j\colon j > i}^n (-b_i)(-b_j)\operatorname{cov}(X_i,X_j) $$ which would be weird, because the predictions are clearly correlated, but I don't see it accounting for any additional covariance. Perhaps it is already accounted for in the individual profile covariances?
Is this accurate? Am I missing something? A proof of this (or disproof) would be welcome.
