# Calculation of IV Standard Errors using R

One thing that I cannot find very clear, perhaps because I am still new in using R, is the derivation of the standard errors in instrumental variable packages using R (e.g., ivreg, sem, plm). I am trying to verify that the standard errors do indeed meet the following requirement found in (Wooldridge, 2009, 511):

$$SE_{IV} = \sqrt{\frac{RSS}{n-K}/\sqrt{TSS_X \times R^2_{X,Z}}}$$

where $X$ is an endogenous regressor and $Z$ an instrumental variable for $X$. The derivation of standard errors for OLS are provided as a convenience below:

$$SE_{OLS} = \sqrt{\frac{RSS}{n-K}/\sqrt{TSS_X}}$$

I've looked at the r code and due to my inability I have not been able to verify standard errors are calculated this way. I've looked at other posts in cross validated and it is not clear from these discussions either. If anyone has some understanding in this area I would be grateful for your input.

When implementing ivreg() I've followed the exposition in Greene (2003) and also compared the results against Stata for which there is an FAQ about 2SLS with some useful details at http://www.stata.com/support/faqs/statistics/two-stage-least-squares/.
Specifically, if $y$ is the response, $X$ the regressor matrix (with both exogenous and endogenous regressors), and $Z$ the instrument matrix (with both exogenous regressors and actual instrument variables), the $P_Z = Z (Z^\top Z)^{-1} Z^\top$. Then, $\hat \beta = (X^\top P_Z X)^{-1} X^\top P_Z y$ and $\text{Cov}(\hat \beta) = (X^\top P_Z X)^{-1} X^\top P_Z \Omega P_Z X (X^\top P_Z X)^{-1}$ where $\Omega = \text{Cov}(y)$ which by default is assumed to be spherical $\Omega = \sigma^2 I$.
Note that the actual computations within ivreg.fit() look different because instead of the matrix operations above, more efficient tools for solving the estimating equations are used. The two stages of least squares are applied through lm.fit() (or lm.wfit() respectively).