st: how many d.f. in the vcv for the within estimator?

["Impavido, Gregorio" <GImpavido@imf.org>]

Dear all, 

I am trying to replicate the result of -xtreg, fe- in mata. Given that the individual effects are not estimated after the within transform, shouldn't the degrees of freedom used for the estimate of the variance of the residuals be (N*T-k-1) instead of N*(T-1)-k (where N=number of panels, T=periods,k=number of regressors)?  I have added the intercept = the meant of the dep variable as in STATA manual.

I paste the code used below which reproduces the results after if force df_r= N*(T-1)-k.  Any suggestion would be welcome.

With kind regards
Gregorio
==================
* this do file reproduces the results of -xtreg, fe- using mata
* it works only with balanced panels
* For corrections and suggestions, Gregorio Impavido (gimpavido@imf.org)
*******************************************************************************
************************************ START ************************************
*******************************************************************************
use "http://www.stata-press.com/data/r9/grunfeld.dta ", clear
rename invest I
rename mvalue F
rename kstock C
sort company time
gen touse=(I!=. & F!=. & C!=.)      // ignore eventual missing obs

* example of panel within estimator
mata:
mata clear                          // clear the workspace

T = 20                              // Number of observations per groups 
N = 10                              // Number of groups 

Z = st_data(.,("F","C"),"touse")    // (NTxk) matrix of regressors
Y = st_data(.,("I"),"touse")        // (NTx1) vector of dep var
i = J(rows(Z),1,1)                  // (NTx1) vector of ones, declare a  
X = Z,i
it = J(T,1,1)                       // (Tx1) vector of ones, declare a 
in = J(N,1,1)                       // (Tx1) vector of ones, declare a  
B = pinv(T)*I(N)#(it*it')           // (NTxNT) between-individual operator
Bbar = pinv(N)*(in*in')#I(T)        // (NTxNT) between-individual operator
Q = I(N*T)-B                        // (NTxNT) Q within transform
* Qbar = I(N*T)-Bbar                  // (NTxNT) Qbar within period operator
Ybar = B*Bbar*Y                     // (NTx1) vector of mean dep var
Zbar = B*Bbar*Z                     // (NTxk) matrix of mean regressor var
Ytilda = Q*Y + Ybar                 // added overall mean
Ztilda = Q*Z + Zbar                 // added overall mean
Xtilda = Ztilda,i                   // (NTxk+1) add column of 1s (intercept)
b_fe = pinv(Xtilda'Xtilda)*Xtilda'*Ytilda	// (k+1x1) vector of beta hat
k = cols(Z)                         // (1x1) No of regressors
u = (Ytilda-Xtilda*b_fe)            // (NTx1) uhat, fitted residuals
df_r = (N*(T-1)-k)                  // (1x1) residual d.f. (shouldn't be NT-k-1??)
rss = (u'*u)                        // (1x1) unrestricted residual sum of squares
mse = rss/df_r                      // (1x1) mean squared error
vcv = mse*pinv(Xtilda'Xtilda)       // (NTxk+1) VCOV matrix
se = sqrt(diagonal(vcv))            // (k+1x1) vector of s.e. of the beta hat
t = b_fe:/se                        // (k+1x1) vector of t statistics
pt = 2*ttail(df_r,abs(t))           // (k+1x1) vector of pvalues
crit = invttail(df_r,0.025)         // (k+1x1) bhat~T(df_r)(b,V(b))
cil = b_fe-crit*se                  // (k+1x1) vector of low CI
cih = b_fe+crit*se                  // (k+1x1) vector of high CI

rss, df_r, mse
b_fe, se, t, pt, cil, cih

end
xtset company time
xtreg I F C, fe                     // to cross check
*******************************************************************************
************************************* END *************************************
*******************************************************************************

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/