Average of Monte-Carlo estimates: Increasing observations vs. iterations

Question

In the context of Monte Carlo simulations, I would like to understand better the difference between increasing the number of iterations vs. the number of observations.

As an example, please consider the following simple R code that is supposed to confirm the unbiasedness of OLS:

set.seed(1)

beta <- c(1, 2)

Niter <- 100
NObs  <- 1000

estimates <- matrix(as.numeric(NA), nrow = Niter, ncol = 2)

for(i in 1:Niter) {

  x <- rnorm(n = NObs)
  u <- rnorm(n = NObs)

  y <- beta[1] + beta[2] * x + u

  reg <- lm(y ~ x)
  estimates[i, ] <- coef(reg)
}

print(colMeans(estimates))
# [1] 1.001352 1.997076

The average of the estimates is pretty close to the true beta values and would become even closer if I increased Niter or NObs.

In the current case, I am only interested in the average of the estimates. If I was interested in their distribution, I would need a sufficiently large Niter to get as many observed estimates as possile. But when I only want to check the average, what is the difference between increasing Niter vs. NObs?

Increasing Niter (100) by 1 gives NObs (1000) observations more. Is this equivalent to increasing NObs by 10 which also gives 1000 additional observations over the 100 iterations?

Answer 1

The expected value of colMeans will allways be unbiased, no matter whether you inrease NObs or Niter, but the difficult part of your question is, whether the variances will change.

We know that estimates $\sim N_2(\mu = \beta, ~ \Sigma = \sigma_\varepsilon^2 \cdot (X_{NObs}^T X_{NObs})^{-1})$ with $\beta$ being beta, $X_{NObs}$ being cbind(1, x) with length(x) == NObs and $\sigma_\varepsilon^2$ being 1. So colMeans(estimates)) is a linear combination of normally distributed random variables. In formula colMeans(estimates)) $\sim N_2(\mu = \beta,~\Sigma = \sigma_\varepsilon^2 \cdot (X_{NObs}^T X_{NObs})^{-1} /Niter )$

In your case $(X_{NObs}^T X_{NObs})^{-1} = (n \cdot \sum(x^2) - (\sum(x))^2)^{-1} \cdot \begin{pmatrix} \sum (x^2) & -\sum(x) \\ -\sum(x) & n \end{pmatrix}$ what can be written as $(n^2 \cdot var(x))^{-1} \cdot \begin{pmatrix} \sum (x^2) & -\sum(x) \\ -\sum(x) & n \end{pmatrix}$ with $var(x)$ being the uncorrected empirical variance - here $1$. For reasons of clarity I use $n = NObs$.

So your question is whether the change in $Niter$ from $100$ to $101$ changes $\begin{pmatrix} \sum (x^2)/n^2 & -\sum(x)/n^2 \\ -\sum(x)/n^2 & 1/n \end{pmatrix}/Niter$ in the same way as a change in $NObs$ from $1000$ to $1010$. If we assume, that $\sum_1^{1010}(x^2) = \sum_1^{1000}(x^2)+ \sum_{1001}^{1010}(x^2) = 1.01 \cdot \sum_{1}^{1000} (x^2)$ and $-\sum_1^{1010}(x) = -\sum_1^{1000}(x) - \sum_{1001}^{1010}(x) = - 1.01 \cdot \sum_{1}^{1000} (x)$, incresing NObs by 10 is equivalent to multiply the covariance matrix by a factor $1.01^{-1}$ (note that $n_{new} = 1010 = 1.01 \cdot n$ for $n=1000$).

So $\begin{pmatrix} \sum (x^2)/n^2 & -\sum(x)/n^2 \\ -\sum(x)/n^2 & 1/n \end{pmatrix}/(100+1)$ is equivalent to $\begin{pmatrix} \sum (x^2)/n^2 & -\sum(x)/n^2 \\ -\sum(x)/n^2 & 1/n \end{pmatrix}\cdot (1.01)^{-1} /100$ (because $1/(100+1) = 1/101 = 1/(1.01 \cdot 100)$).

---

I did a small "meta"-simulation study, where I ran your simulation 100 times under both Settings. Setting 1 was to increase Niter by 1, and setting 2 was to increase NObs by 10. In each meta-Simulation run, I saved the colMeans in a matrix colmeans1 and colmeans2. Plotting `the entries from both matrices doesn't contradict the theoretical finding of both settings being equivalent.

set.seed(1)
S <- 100

colmeans1  <- matrix(NA, nrow = S, ncol = 2)
colmeans2  <- matrix(NA, nrow = S, ncol = 2)

#Setting 1: increasing Niter by 1
for(s in 1:S){
beta <- c(1, 2)

Niter <- 100 + 1
NObs  <- 1000

estimates <- matrix(as.numeric(NA), nrow = Niter, ncol = 2)

for(i in 1:Niter) {

  x <- rnorm(n = NObs)
  u <- rnorm(n = NObs)

  y <- beta[1] + beta[2] * x + u

  reg <- lm(y ~ x)
  estimates[i, ] <- coef(reg)
}
colmeans1[s,]  <- colMeans(estimates)
}



#Setting 1: increasing NObs by 10
for(s in 1:S){
beta <- c(1, 2)

Niter <- 100 
NObs  <- 1000 + 10

estimates <- matrix(as.numeric(NA), nrow = Niter, ncol = 2)

for(i in 1:Niter) {

  x <- rnorm(n = NObs)
  u <- rnorm(n = NObs)

  y <- beta[1] + beta[2] * x + u

  reg <- lm(y ~ x)
  estimates[i, ] <- coef(reg)
}
colmeans2[s,]  <- colMeans(estimates)
}