How can I code this equation with double summation in R?

Question

So I'm having hard time coding the above equation, mainly the part which contains that double sum over i's and over j.

I'n my case, my n = 200 and p = 15. My yi:s are in a vector Y = (y1,y2,...yn) that is vector of length 200 and Xij:s are in a matrix which has 15 columns and 200 rows. Bj:s are in a vector of length 15.

My own solution, which I'm fairly certain is wrong, is this:

b0 <- 1/200 * sum(Y - sum(matr*b))

And here is code which you can use to reproduce my vectors and matrix:

matr <- t(mvrnorm(15,mu= rep(0,200),diag(1,nrow = 200)))
Y <- rnorm(n = 200)
b <- rnorm(n = 15)

Answer 1

Use matrix multiplication:

mean(y - x %*% b)

Note that if y and x are known and b is the least squares regression estimate of the coefficients then we can write it as:

fm <- lm(y ~ x + 0)
mean(resid(fm))

and that necessarily equals 0 if there is an intercept, i.e. a constant column in x, since the residual vector must be orthogonal to the range of x and taking the mean is the same as taking the inner product of the residuals and a vector whose elements are all the same (and equal to 1/n).