How to use a for loop in this Bernoulli exercise in R? I have a (for me) very difficult exercise, where I want to ask you guys for help as I am still a newbie on R. (But I'm getting there!)
So my question looks like this:
Given Bernoulli distribution with p=0.78
For each sample size n simulate r = 10,000 draws with a for loop from that Bernoulli distribution with p = 0.78 then calculate the standardized sample mean for each of these draw. Replicate all four panels from Figure 2.9 See Figure 2.9 attached to this.
So I know I have to use
rbinom()
for this.
But that's it. I really appreciate you guys help on this one.

(From Stock & Watson (2015), Introduction to Econometrics, 3rd Ed., Fig. 2.9)
 A: The instruction in the exercise to use a loop is bad advice.  The rbinom function is already capable of simulating vectors of values, so there is no need for a loop.  The simplest thing to do here is to create an $r \times N$ matrix of simulated Bernoulli random variables taking $N=100$ so that you have enough sample size to meet the requirements of the four specified values of $n$.  Assuming you don't mind having nested simulations (which is not a problem) you then have all the simulated values you need to construct the figures.  Here is some simple code to generate a "reproducible" matrix of simulated Bernoulli random variables.  (Note that you can also use the sample.int function to efficiently simulate Bernoulli random variables.)  It will not give you the exact outcomes used in the graphs in that book, but it will nonetheless give you reproducible simulated values.
#Set parameters
N    <- 100
r    <- 10000
PROB <- 0.78

#Create matrix for simulated values
set.seed(1)
SIMULATIONS <- matrix(rbinom(r*N, size = 1, prob = PROB), nrow = r, ncol = N)
colnames(SIMULATIONS) <- sprintf('Sample[%s]', 1:N)

Once you have simulated the matrix SIMULATIONS you then have $r$ rows of simulated values for $N$ sample values.  You can obtain the relevant simulations of the sample means by using the rowMeans function on the relevant subsets of the matrix.  You can then use appropriate plotting functions to construct the required plots.  This gives similar results to the graphs you have shown.
#Create matrix of standardised sample means
STD.MEANS <- matrix(0, nrow = R, ncol = 4)
colnames(MEANS) <- c('n[2]', 'n[5]', 'n[25]', 'n[100]')
STD.MEANS[, 1] <-   sqrt(2)*(rowMeans(SIMULATIONS[, 1:2])   - PROB)/sqrt(PROB*(1-PROB))
STD.MEANS[, 2] <-   sqrt(5)*(rowMeans(SIMULATIONS[, 1:5])   - PROB)/sqrt(PROB*(1-PROB))
STD.MEANS[, 3] <-  sqrt(25)*(rowMeans(SIMULATIONS[, 1:25])  - PROB)/sqrt(PROB*(1-PROB))
STD.MEANS[, 4] <- sqrt(100)*(rowMeans(SIMULATIONS[, 1:100]) - PROB)/sqrt(PROB*(1-PROB))

#Plot the histograms
par(mfrow = c(2,2))
hist(STD.MEANS[, 1], prob = TRUE, col = "skyblue2", xlim = c(-5, 5), 
     main = '(n = 2)', xlab = 'Standardised Sample Mean')
curve(dnorm(x), add = TRUE, lwd = 2)
hist(STD.MEANS[, 2], prob = TRUE, col = "skyblue2", xlim = c(-5, 5), 
     main = '(n = 5)', xlab = 'Standardised Sample Mean')
curve(dnorm(x), add = TRUE, lwd = 2)
hist(STD.MEANS[, 3], prob = TRUE, col = "skyblue2", xlim = c(-5, 5), 
     main = '(n = 25)', xlab = 'Standardised Sample Mean')
curve(dnorm(x), add = TRUE, lwd = 2)
hist(STD.MEANS[, 4], prob = TRUE, col = "skyblue2", xlim = c(-5, 5), 
     main = '(n = 100)', xlab = 'Standardised Sample Mean')
curve(dnorm(x), add = TRUE, lwd = 2)


A: To begin, I agree with @Ben's(+1) statements about avoiding explicit
loops when possible. I have used for loops because they seem
to be required for your exercise.
Standardization is done outside the for loop, using means and standard deviations from the 10,000 averages a.
Here is a simulation in R for the case $n = 25.$
set.seed(121)
n = 25;  p = 0.78
r = 10^4;  a = numeric(r)
for(i in 1:r) {
  a[i] = mean(rbinom(n, 1, .78))  
  }
mean(a);  sd(a)
z = (a-mean(a))/sd(a)
cp = seq(-5.75, 5.75, length=13)
hdr = "n=25: Standardized Value of Sample Average"
hist(z, prob=T, br=cp, ylim=c(0,.4), col="skyblue2", main=hdr)
 curve(dnorm(x), add=T, lwd=2)


Note about making histograms: I have used standard graphics from the base of R to make
the histogram and superimpose the standard normal density curve.
Even though $r = 10,000$ means have been generated so that a
has $r$ values, there are not many unique values in a--fifteen
in my simulation (some relatively rare). If the tied values are not equitably proportioned
among the histogram bins, you get some strange-looking histograms
that make the values of z look far from normal. By choosing 13 bins
I got a nice plot. (Roughly speaking, there are two z-values per bin, with some empty bins at the ends.)
length(unique(a))
[1] 15
table(a)
a
0.44 0.48 0.52 0.56  0.6 0.64 0.68 0.72 0.76  0.8 0.84 0.88 0.92 0.96    1 
   1   12   22   86  212  459  847 1381 1795 1914 1591 1037  495  123   25 

