Is bootstrapping appropriate for estimating a multivariate normal covariance matrix using a small sample size? Let $p$ be a matrix, where each row represents an observation of a 4-variate normally distributed random variable $\mathcal{N}_4(\mu,\Sigma)$. 


*

*Is there any Bootstrap methode to get a good estimation for Σ?

*If not, is the number of the following sample enough to bootstrap the distribution of any statistic T which operates on the population where Y comes from? 
Here's what I tried so far:
Y<-data.frame(response=c(10,19,27,28,9,13,25,29,4,10,20,18,5,6,12,17),
               treatment=factor(rep(1:4,4)),
               subject=factor(rep(1:4,each=4))
               )

p<-matrix(Y$response,4,4,byrow=T)
B<-1000
sampleB<-sample(1:4,4*B,replace=T)
fit<-lm(p[sampleB,]~1)
cov(residuals(fit))

I also tried 
require(nlme)
require(mgcv)
nSubj <- 20
sampleB<-sample(1:4,nSubj,replace=T)
y<-data.frame(response=c(t(p[sampleB,])),
           treatment=factor(rep(1:4,nSubj)),
           subject=factor(rep(1:nSubj,each=4))
           )

fit <- lme(response~-1+treatment,y,random=~1|subject,correlation=corSymm())
extract.lme.cov(fit,y)[1:4,1:4]

but I get the error code: 
Error in lme.formula(response ~ -1 + treatment, y, random = ~1 | subject,  : 
nlminb problem, convergence error code = 1
message = iteration limit reached without convergence (10)

 A: 
Is the Bootstrap a good option to estimate $\Sigma$? 

No, the boostrap will help you infer about the uncertainty of your sample estimate. Specifically, it might be used to get confidence intervals on the elements of $\widehat{\Sigma}$.

And how does the bootstrap work?

The approach is to create $R$ replicate datasets from the original dataset by resampling the observations with replacement. Then you compute the estimate of interest on each of the $R$ replicates, in your case the covariance matrix, for each of the $R$ replicates, obtaining $\widehat{\Sigma}^1, \ldots, \widehat{\Sigma}^R$. Confidence intervals for $\widehat{\Sigma}$ can then be computed empiricaly from $\widehat{\Sigma}^1, \ldots, \widehat{\Sigma}^R$.
For further information, you might want to have a look at the Wikipedia page.

Edit

Could I use $\widetilde{\Sigma}_R = R^{-1} \sum_{r=1}^R \widehat{\Sigma}^r$ instead of $\widehat{\Sigma}$?

Actually, the matrix $\widetilde{\Sigma}_R$ is an estimate of ${\rm E} (\widehat{\Sigma}^r)$, where $\widehat{\Sigma}^r$ is the estimate of the covariance matrix based on a bootstrap replication of the initial sample. 
Bootstrapping comes down to sample from the empirical distribution $\widehat{F}$.
Therefore, I think, but I don't have a formal proof, that $\widetilde{\Sigma}_R = {\rm E} (\widehat{\Sigma}^r) \to \widehat{\Sigma}$ as $R \to \infty$.
So, I think you could use $\widetilde{\Sigma}_R$ instead of $\widehat{\Sigma}$, but that would be like using a sledgehammer to crack a nut.
The R code below is a numerical investigation, which is by no means of proof of the above assertion about the convergence. The structure of dependence is a Gumbel copula, and the margins are two standard normal distribution.
## Initialization
library(copula)
set.seed(531)
n <- 200            # Number of observations in the original sample
R <- 10000          # Number of replications
## Specification for the dependence structure (Gumbel copula)
spec.cl <- archmCopula("gumbel", 1.2)
## Create a fake original dataset
pseudo  <- rCopula(n, spec.cl)
obs     <- qnorm(pseudo)
cov.obs <- cov(obs)[1, 2]
## Get an idea of the "true" covariance
pseudo   <- rCopula(10000, spec.cl)
obs.big  <- qnorm(pseudo)
cov.true <- cov(obs.big)[1, 2]
## Get the bootstrap covariances
cov.sim <- sapply(1:R,
                  function(i, x, n){x.boot <- x[sample(1:n, size = n, replace = TRUE), ]
                                    cov(x.boot)[1, 2]},
                  x = obs, n = n)
## Visualization
plot(1:R, cov.sim, xlab = "Replication", ylab = "", pch = 16, cex = 0.7, col ="grey",
     ylim = quantile(cov.sim, probs = c(0.1, 0.9)))
lines(1:R, rep(cov.true, R), col = "green", lwd = 2)
lines(1:R, rep(cov.obs, R), col = "red", lwd = 2)
lines(1:R, cumsum(cov.sim)/(1:R), col = "blue", lwd = 2)
legend("topright", legend = c("Boot cov", "True", "Initial", "Boot average"),
       col = c("grey", "green", "red", "blue"),
       bg = "white", pch = c(16, NA, NA, NA), lwd = c(NA, 2, 2, 2))

The blue line corresponds to $\widetilde{\Sigma}_R$ as a function of the number of replications, and the red line is $\widehat{\Sigma}$.

