Underestimation of standard error

Question

The formula to calculate standard error (aka standard deviation of the sample mean) is $\sigma_{\bar{x}} = \frac{ \sigma} {\sqrt{n}}$. But we generally don't know the population standard deviation $\sigma$. So we have to first estimate $\sigma$ from our sample (note: it is a biased estimator). Using this estimated standard deviation, our new formula is $\sigma_{\bar{x}} \approx \frac{ \sigma_x} {\sqrt{n}}$.

I did a Monte Carlo simulation, and found that this estimated standard error is systematically lower than the exact standard error. I am wondering why this is the case? It will be much appreciated if anyone can prove this.

Here I attached some R code if you want to try. Notice that the 3rd result is lower than the others.

# draw 100000 samples, each sample has 5 data point
mat = matrix(nrow = 100000,ncol = 5)
for (i in 1:100000){
  mat[i,1:5] <- rnorm(5)
}



# Method1: Calculate standard error by definition: standard error is the standard deviation of sample means
sd(rowMeans(mat))


# Method2: Calculate standard error using the population std: sigma/sqrt(n)
1/sqrt(5)


# Method3: Calculate standard error using the sample std: sigma_hat/sqrt(n)
row_std = apply(mat, 1, sd)
row_se = row_std/sqrt(5)
mean(row_se)

Answer 1

This is what's going on

> summary(row_std)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.05212 0.69485 0.91762 0.94109 1.15915 2.56883 
> summary(row_std^2)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.002717 0.482820 0.842034 1.001564 1.343628 6.598870 

Because $s^2$ is unbiased for $\sigma^2$, it is not possible for $s$ to be unbiased for $\sigma$.

To a first approximation $$E[f(Z)]=f(E[Z])+f''(E[Z])\dfrac{\mathrm{var}(Z)}{2}$$ where the second derivative is $-\frac{1}{2}\sigma^{-3}$ and the variance is $\frac{\sigma^2}{n}\times (2-\frac{2}{n-1})$ (for a Gaussian), giving a first correction of -0.075. This overcorrects. Also, it depends on the unknown $\sigma^2$ and kurtosis, so it wouldn't be estimated all that well from five observations if we weren't pretending we knew the values.

More importantly, if the data are Gaussian, the distribution of $s^2$ is already taken into account in computing the confidence interval

> in_interval<-function(theta,hat,se,tcrit){ (hat-tcrit*se <= theta) & (hat+tcrit*se>=theta)}
> meanhat<-rowMeans(mat)
> table(in_interval(0,meanhat, row_se,abs(qt(.025,4))))

FALSE  TRUE 
 4783 95217 

Answer 2

By Jensen’s inequality, the sample standard deviation is an underestimate (in expectation) of the true standard deviation, since the square root is concave and $S^2$ is unbiased for the second central moment.