You can even use outer
and diag
to get the same result
> p <- drop(p_hat)
> variance <- p*(1-p)
> covariance <- -outer(p, p)
> diag(covariance) <- variance
> covariance
[,1] [,2] [,3]
[1,] 0.090 -0.0290 -0.0610
[2,] -0.029 0.2059 -0.1769
[3,] -0.061 -0.1769 0.2379