Estimated variance using linear factor models This might sound like a trivial question, but I just can't get my head around it! 
I am trying to evaluate the performance of different risk models (linear factor models), by computing the bias statistic. It is simply a function of estimated variance of the asset you are testing for and its forward return.
Let's consider the following model, $M_1$
$R$ = a + $b_1$$F_1$ + $b_2$$F_2$ + e 
= a + $B*F$ + e
where R is the return of the asset, F are the factor premiums, B are the factor exposures, a is a constant and e is the residual
The estimated asset variance using $M_1$ is:
$var(R)$ = $B^T$$\Sigma$$B$ + $\Delta$
where $\Sigma$ is the variance-covariance matrix of factors F, and $\Delta$ is the variance of the regression residuals.
However, if I consider model $M_2$, which has completely different factors that $M_1$
$R$ = a + $b_3$$F_3$ + $b_4$$F_4$ + e 
and compute the estimated variance of the asset, I get the exact same value!!
Here is the link to a snapshot of the data: data. The first column (after the date column) is the dependent variable, where the rest of the columns are factors $F_1$, $F_2$, $F_3$, $F_4$. Here is the code I used to compute the forecasted risk:
library(zoo)

data = read.csv(csv, stringsAsFactors = FALSE, na.strings= "NA")
data1.zoo = zoo(data[,c(2,3,4)], order.by = as.Date(data[,1], format="%Y-%m-%d"))
data2.zoo = zoo(data[,c(2,5,6)], order.by = as.Date(data[,1], format="%Y-%m-%d"))

Function to estimate forecasted variance
get_est_var <- function(z) {
   df = as.data.frame(z)
   reg = lm(sp500 ~., data = as.data.frame(z), na.action = na.omit)
   beta = t(as.matrix(reg$coefficients[-1]))
   res = as.matrix(reg$residuals)
   sigma = as.matrix(cov(df[,-1]))

   est.var = (beta %*% tcrossprod(sigma, beta)) + (deviance(reg)/df.residual(reg))   
   est.var
}

risk1 = sqrt(get_est_var(data1.zoo))
risk2 = sqrt(get_est_var(data2.zoo))

You will see that the values risk1 and risk2 are pretty similar although the factors used are completely different!!!
Can anyone please enlighten me as to how do we compute model specific estimated asset variance? Thanks for your help!
REFERENCE
Here is a reference regarding the computation of variance in the context of factor models: link. Please take a look at the formula at the bottom of page 7. The covariance of the factors formula is on top of page 4. Please note that every authoritative source on portfolio management corroborates these formulas.
Here is another reference about the implementation of risk models in R: link to implementation in R. Please see page 10 and page 12.
 A: There exists a basic result in least squares algebra that says that, in a linear regression and using OLS estimation, regressing a variable on a set of regressors plus a constant, is equivalent (in a mathematically exact way) to estimating using the dependent variable and the regressors as deviations from their respective means.
The specification that is presented in the reference you provided, adjusted for univariate  time-series regression with a sample of size $T$ (I am sorry but I always use $T$ to denote the time-dimension, unlike the reference), with one dependent variable is (in vector matrix notation to represent the whole sample)
$$\mathbf R = a \cdot\mathbf 1 + \mathbf F\mathbf B + \mathbf u$$
Now write $\mathbf D = I_{T\times T} -\mathbf 1\left(\mathbf 1'\mathbf1\right)^{-1}\mathbf1' $ for the special matrix that de-means a column vector, and define
$$\tilde {\mathbf R} = \mathbf D \mathbf R,\;\; \tilde {\mathbf F} = \mathbf D \mathbf F$$
Then you have equivalently
$$\tilde {\mathbf R} = \tilde {\mathbf F}\mathbf B + \mathbf u$$
From this alternative (but again, equivalent) representation of the model we can see that the sample variance of the dependent variable (spare me the bias correction term) can be expressed as (using estimated magnitudes)
$$\operatorname{\hat Var}(R) = \frac 1n \tilde {\mathbf R}'\tilde {\mathbf R} = \frac 1n\left[\tilde {\mathbf F}\hat {\mathbf B} + \hat {\mathbf u}\right]'\left[\tilde {\mathbf F}\hat {\mathbf B} + \hat {\mathbf u}\right]$$
$$=\frac 1n\hat {\mathbf B}'\tilde {\mathbf F}'\tilde {\mathbf F}\hat {\mathbf B} + \frac 1n\hat {\mathbf B}'\tilde {\mathbf F}'\hat {\mathbf u} + \frac 1n\hat {\mathbf u}'\tilde {\mathbf F}\hat {\mathbf B}+\frac 1n\hat {\mathbf u}'\hat {\mathbf u}$$
Note the following: by construction, the regressors are orthogonal to the error term, and hence the 2nd and 3d terms are exactly zero. Moreover, $\frac 1n\tilde {\mathbf F}'\tilde {\mathbf F} = \Sigma$ (the $\Sigma$ in the question) so 
$$\operatorname{\hat Var}(R) = \hat {\mathbf B}'\Sigma\hat {\mathbf B} + \frac 1n\hat {\mathbf u}'\hat {\mathbf u}\qquad [1]$$
which is the expression for the sample variance of $R$ given in the reference, which we see that it can be obtained easily through this alternative representation of the model.
The LHS of this equation is the sample variance of the dependent variable. This magnitude should not be affected by the choice of regressors -and so taking any set of regressors, expressing them in mean deviation form, running a regression, obtaining a new set of residuals and then plugging all these into the RHS of the equation, should give exactly the same result. And it does, as the OP found numerically, which can also be shown to hold for any set of regressors. Define $\mathbf M_F$ to be the residual-maker matrix of the regression, $\mathbf M_F = I - \mathbf P_F$ where $\mathbf P_F$ is the projection matrix $P_F = \tilde {\mathbf F}\Big(\tilde {\mathbf F}'\tilde {\mathbf F}\Big)^{-1}\tilde {\mathbf F}'$
Then, using the expressions for the estimated coefficients and the residuals we have
$$n\cdot RHS =\left[\Big(\tilde {\mathbf F}'\tilde {\mathbf F}\Big)^{-1}\tilde {\mathbf F}'\tilde {\mathbf R}\right]'  \tilde {\mathbf F}'\tilde {\mathbf F}\Big(\tilde {\mathbf F}'\tilde {\mathbf F}\Big)^{-1}\tilde {\mathbf F}'\tilde {\mathbf R} + \left(\mathbf M_F\tilde {\mathbf R}\right)'\left(\mathbf M_F\tilde {\mathbf R}\right)$$
$$= \tilde {\mathbf R}'\tilde {\mathbf F}\Big(\tilde {\mathbf F}'\tilde {\mathbf F}\Big)^{-1}\tilde {\mathbf F}'\tilde {\mathbf F}\Big(\tilde {\mathbf F}'\tilde {\mathbf F}\Big)^{-1}\tilde {\mathbf F}'\tilde {\mathbf R} + \tilde {\mathbf R}'\mathbf M_F\mathbf M_F\tilde {\mathbf R}$$
where we have used the fact that $\mathbf M_F$ is always a symmetric matrix. Simplifying and using the fact that the residual maker matrix is also idempotent, we have
$$n\cdot RHS =\tilde {\mathbf R}'\tilde {\mathbf F}\Big(\tilde {\mathbf F}'\tilde {\mathbf F}\Big)^{-1}\tilde {\mathbf F}'\tilde {\mathbf R} + \tilde {\mathbf R}'\mathbf M_F\tilde {\mathbf R}$$
and using the relation between the residula maker and the projection matrix we have
$$n\cdot RHS =\tilde {\mathbf R}'\mathbf P_F\tilde {\mathbf R} + \tilde {\mathbf R}'\mathbf M_F\tilde {\mathbf R} = \tilde {\mathbf R}'(I-\mathbf M_F)\tilde {\mathbf R} + \tilde {\mathbf R}'\mathbf M_F\tilde {\mathbf R}$$
$$\Rightarrow n\cdot RHS = \tilde {\mathbf R}'\tilde {\mathbf R} $$
$$\Rightarrow RHS = \frac 1n \tilde {\mathbf R}'\tilde {\mathbf R} $$
So the RHS of eq. $[1]$ is composed in such a way as to be mathematically equal with the result we would obtain if we simply calculated the sample variance of the dependent variable, irrespective of the regressors chosen.  
What the choice of the regressors affects is the allocation of the sample variance into a "common factor component" $\hat {\mathbf B}'\Sigma\hat {\mathbf B}$ and into an "asset specific" component $\frac 1n\hat {\mathbf u}'\hat {\mathbf u}$ (which is the "translation" in the context of the specific model of the traditional statement about the"explained" and "unexplained" portion of the variance of the dependent variable).  
So, you want to compute something called the "bias statistic", that involves "the estimated variance of the asset and its forward return". If you should  expect that this statistic should have a different value depending on the regressors chosen, then, in light of the above, this could happen if
1) In the bias statistic, only one of the two parts of the decomposed sample variance enters(i.e. only one of the two terms of the RHS of eq. $[1]$  
and / or
2) The choice of regressors affects the "forward return". 
I am not familiar with the definition of the bias statistic and what it attempts to measure so I cannot help you further than that.
A: The values are not exactly the same, but are very close.  This appears to be purely accidental when we look at what goes into them:
Model    beta'*Sigma*beta  Deviance     DF    Deviance/DF   "est.var"
    1    0.0001504468      0.0005199823 1256  4.139987e-07  0.0001508608
    2    9.168255e-08      0.1896668    1256  0.0001510086  0.0001511003

In the first model the result is approximately $0.000151$ because $\beta'\Sigma\beta$ is almost this value and the deviance (divided by the DF) adds essentially nothing.  In the second model $\beta'\Sigma\beta$ contributes almost nothing but the deviance (divided by the DF) is approximately $0.000151.$
It is a very good idea to look into results that seem like more than coincidence, as you have done here, but accidents do happen.
