Standard error of regression coefficient without raw data After searching here:
Perform simple regression without raw data
I am still curious about this.
Is it possible to derive the standard error of a regression coefficient from summary data alone? 
E.g., assume we are given the following variance-covariance matrix.
$\begin{bmatrix}Var(X) & Cov(X,Y)\\Cov(X,Y) & Var(Y)\end{bmatrix}$
We can derive the regression coefficient $\beta_{XY} = Cov(X,Y) / Var(X)$.
Given a specific n, is it possible to derive the standard error of $\beta$ as well? If so, which formula is being used? It appears that all formulas for regression standard errors that I could find assume that you know the variance of residuals of the regression, which we don't know from summary data alone.
---UPDATE---
Based on feedback by Greg and whuber, I am now at a point, where I know that
$\sigma^2(b) = \sigma^2(X'X)^{-1}=\sigma^2(Var(X)(n-1))^{-1}$,
but typically matrix formulas for the standard error of $b$ then note that $\sigma^2$, the unknown variance of the errors, is estimated with MSE. 
Based on comments by whuber, sum of squares residuals is:
$Y'Y - bX'X$, which according to Greg can be spelled out as
$Var(Y)(n-1) - b(Var(x)(n-1))$
And to go from here to MSE, I think (?) all that is left is divide by n-2,
so
$(Var(Y)(n-1) - b(Var(x)(n-1))) / (n-2)$
And plugging this all in yields:
$\sigma^2(b) =((Var(Y)(n-1) - b(Var(X)(n-1))) / (n-2))\times(Var(X)(n-1))^{-1}$
I will try to see if this works out in R, and will report back.
If anyone sees something blatantly wrong, I appreciate the heads up.
 A: One quick and dirty approach when you have the summaries but not the raw data is to generate a dataset with the specified summaries and sample size and then run the simulated data through your regular routines to compute whatever you want.  The mvrnorm function in the MASS package for R will generate random normal data with a given mean vector and covariance matrix.  I am sure other programs will as well (or you can create your own by generating data and then multiplying by the appropriate matrix).
From the theoretical side you can start by thinking about having all your variables centered (mean subtracted so that the current mean is exactly 0, this really only affects the intercept) then $x'x=var(x)\times (n-1)$, $x'y=cov(x,y)\times(n-1)$, and $y'y=var(y)\times(n-1)$.  So just plug those values into the matrix version of the formula for the standard error that you want.
A: With the usual notation, arrange the independent variables by columns in an $n\times (p+1)$ array $X$, with one of them filled with the constant $1$, and arrange the dependent values $Y$ into an $n$-vector (also a column).  Assuming the appropriate inverses exist, the desired formulas are
$$b = (X^\prime X)^{-1} X^\prime Y$$
for the $p+1$ parameter estimates $b$,
$$s^2 = (Y^\prime Y - b^\prime X^\prime Y)/(n-p-1)$$
for the residual standard error, and
$$V = (X^\prime X)^{-1} s^2$$
for the variance-covariance matrix of $b$.
When only the summary statistics $\newcommand{\m}{\mathrm m} \m_X$ (the means of the columns of $X$, forming a $p+1$-covector which includes a mean of $1$ for the constant), $m_Y$ (the mean of $Y$), $\text{Cov}(Y)$, $\text{Var}(Y)$, and $\text{Cov}(X,Y)$ are available, first recover the needed values via
$$(X^\prime X)_0 = (n-1)\text{Cov}(X) + n \m_X \m_X^\prime,$$
$$Y^\prime Y = (n-1)\text{Var}(Y) + n m_Y,$$
$$(X^\prime Y)_0 = (n-1)\text{Cov}(X,Y) + n m_Y \m_X.$$
Then to obtain $X^\prime X$, border $(X^\prime X)_0$ symmetrically by a vector of column sums (given by $n \m_X$) with the value $n$ on the diagonal;  and to obtain $X^\prime Y$, augment the vector $(X^\prime Y)_0$ with the sum $n m_Y$.  For instance, when using the first column for the constant term, these bordered matrices will look like
$$X^\prime X = \pmatrix{
  n & n\m_X \\
  n\m_X^\prime & (X^\prime X)_0
}$$
and
$$X^\prime Y = \left(n m_Y, (X^\prime Y)_0\right)^\prime$$
in block-matrix form.
If the means are not available--the question did not indicate they are--then replace them all with zeros.  The output will estimate an "intercept" of $0$, of course, and its standard error of the intercept will likely be incorrect, but the remaining coefficient estimates and standard errors will be correct.

Code
The following R code generates data, uses the preceding formulas to compute $b$, $s^2$, and the diagonal of $V$ from only the means and covariances of the data (along with the values of $n$ and $p$ of course), and compares them to standard least-squares output derived from the data.  In all examples I have run (including multiple regression with $p\gt 1$) agreement is exact to the default output precision (about seven decimal places).
For simplicity--to avoid doing essentially the same set of operations three times--this code first combines all the summary data into a single matrix v and then extracts $X^\prime X$, $X^\prime Y$, and $Y^\prime Y$ from its entries.  The comments note what is happening at each step.
n <- 24
p <- 3
beta <- seq(-p, p, length.out=p)# The model
set.seed(17)
x <- matrix(rnorm(n*p), ncol=p) # Independent variables
y <- x %*% beta + rnorm(n)      # Dependent variable plus error
#
# Compute the first and second order data summaries.
#
m <- rep(0, p+1)                # Default means
m <- colMeans(cbind(x,y))       # If means are available--comment out otherwise
v <- cov(cbind(x,y))            # All variances and covariances
# 
# From this point on, only the summaries `m` and `v` are used for the calculations
# (along with `n` and `p`, of course).
#
m <- m * n                      # Compute column sums
v <- v * (n-1)                  # Recover sums of squares of residuals
v <- v + outer(m, m)/n          # Adjust to obtain the sums of squares
v <- rbind(c(n, m), cbind(m, v))# Border with the sums and the data count
xx <- v[-(p+2), -(p+2)]         # Extract X'X
xy <- v[-(p+2), p+2]            # Extract X'Y
yy <- v[p+2, p+2]               # Extract Y'Y
b <- solve(xx, xy)              # Compute the coefficient estimates
s2 <- (yy - b %*% xy) / (n-p-1) # Compute the residual variance estimate
#
# Compare to `lm`.
#
fit <- summary(lm(y ~ x))
(rbind(Correct=coef(fit)[, "Estimate"], From.summary=b))    # Coeff. estimates
(c(Correct=fit$sigma, From.summary=sqrt(s2)))               # Residual SE
#
# The SE of the intercept will be incorrect unless true means are provided.
#
se <- sqrt(diag(solve(xx) * c(s2))) # Remove `diag` to compute the full var-covar matrix
(rbind(Correct=coef(fit)[, "Std. Error"], From.summary=se)) # Coeff. SEs

