Decomposition of MSE, an example of Standard Deviation

Given such a basic setting (a setting that appeared in my Research Module of Econometrics and Statistics):

Here sigma_hat is a biased estimator, and sigma_tild is an unbiased estimator of standard deviation, but actually sigma_tild as standard deviation is also biased, according to https://en.wikipedia.org/wiki/Variance, thus here I also show the Bias^2 of sigma_tild^2, which is the variance.

Here as shown above, the blue lines are the Bias^2 of sigma_hat and sigma_tild, and we can see that just like normally believed, the "unbiased" one (which is called Bias^2 (n - 1) ) really has a lower bias square than the "biased" one. However, as the sample size "n_seq" is small, they are all biased. Nevertheless, the green triangles, which stand for the unbiased variance estimator, keeps being really unbiased, which equals zero.

Meanwhile, the MSE of the two estimators are similar, and we can observe that MSE = Var + Bias^2.

Code:

muf <- function(X,n){
  output <- sum(X)/n
  return(output)
}

sigmafn <- function(X,n){
  dev <- X-mean(X)
  dev <- 1/n*sum(dev^2)
  output <- sqrt(dev)
  return(output)
}

sigmafn1 <- function(X,n){
  dev <- X-mean(X)
  dev <- 1/(n-1)*sum(dev^2)
  output <- sqrt(dev)
  return(output)
}

sigma2fn1 <- function(X,n){
  dev <- X-mean(X)
  dev <- 1/(n-1)*sum(dev^2)
  return(dev)
}

n_seq <- seq(2,30)
mu=1
sigma <- 1
Times <- 10000

resultsmu <- matrix(NA,Times,length(n_seq))
resultsigma <- matrix(NA,Times,length(n_seq))
resultsigman <- matrix(NA,Times,length(n_seq))
resultsigman2 <- matrix(NA,Times,length(n_seq))

set.seed(527)
for (i in 1:Times){
  for (j in 1:length(n_seq)) {
  X <- rnorm(n_seq[j],1,1)
  resultsmu[i,j] <- muf(X,n_seq[j])
  resultsigma[i,j] <- sigmafn(X,n_seq[j])
  resultsigman[i,j] <- sigmafn1(X,n_seq[j])
  resultsigman2[i,j] <- sigma2fn1(X,n_seq[j])
  }
}

MSEf <- function(si,Ti){
  dev2 <- sum((si-1)^2)
  output <- dev2/Ti
  return(output)
}

Varf <- function(si,Ti){
  mu <- mean(si)
  dev2 <- sum((si-mu)^2)
  out <- dev2/Ti
  return(out)
}

Bias2f <- function(si,Ti){
  bia <- (mean(si)-1)^2
  return(bia)
}


MSE <- c()
Var <- c()
Bias2 <- c()

MSEn <- c()
Varn <- c()
Bias2n <- c()
Bias22 <- c()

for (i in 1:length(n_seq)) {
  MSE[i] <- MSEf(resultsigma[,i],Times)
  MSEn[i] <- MSEf(resultsigman[,i],Times)
  Var[i] <- Varf(resultsigma[,i],Times)
  Varn[i] <- Varf(resultsigman[,i],Times)
  Bias2[i] <- Bias2f(resultsigma[,i],Times)
  Bias2n[i] <- Bias2f(resultsigman[,i],Times)
  Bias22[i] <- Bias2f(resultsigman2[,i],Times)
}

plot(n_seq,Bias22)
lines(n_seq,Bias2n)


plot(n_seq,MSE,ylim = c(0,0.5),type="l",
     main = "Decomposition of MSE",lwd=2)
lines(n_seq,Var,col="red",lwd=2)
lines(n_seq,Bias2,col="blue",lwd=2)
points(n_seq,MSEn,pch=2,"o")
points(n_seq,Varn,col="red",pch=3,"o")
points(n_seq,Bias2n,col="blue",pch=2,"o")
points(n_seq,Bias22,col="green",pch=6,"o")
abline(0,0,lty=2)
legend("topright",c("MSE","Var","Bias^2",
       "MSE (n-1)","Var (n-1)","Bias^2 (n-1)","Bias^2 of Variance"),
       col=c("black","red","blue","black","red","blue","green"),
       pch=c(19,19,19,2,3,2,6),lwd=c(2,2,2,1,1,1,1))