This kind of problem can be programmed straightforwardly, but that usually results in a tangle of for loops. Some linear algebraic manipulations go a long way in drastically reducing code length (if your computer's memory can handle it).
Consider the original regression equation $$ Y_i= \boldsymbol{X}_i'\boldsymbol{\beta} + \varepsilon_i;\;i=1, \ldots, n $$ where $\boldsymbol{X}_i$ is a $k\times 1$ vector of regressors.You can write this more compactly as $$ \boldsymbol{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon} $$
Then the model that you want to estimate can be stacked up as $$ \begin{bmatrix} \boldsymbol{Y}^{(-1)} \\ \vdots \\ \boldsymbol{Y}^{(-n)} \\ \end{bmatrix} = \begin{bmatrix} \mathbf{X}^{(-1)} & & \\ & \ddots & \\ & & \mathbf{X}^{(-n)} \end{bmatrix} \underbrace{\begin{bmatrix} \boldsymbol{\beta}^{(-1)} \\ \vdots \\ \boldsymbol{\beta}^{(-n)} \\ \end{bmatrix}}_{\equiv \boldsymbol{\beta}^-} + \begin{bmatrix} \boldsymbol{\varepsilon}^{(-1)} \\ \vdots \\ \boldsymbol{\varepsilon}^{(-n)} \\ \end{bmatrix} $$ where $\boldsymbol{Y}^{(-i)} = [Y_1, \ldots, Y_{i-1}, Y_{i+1}, \ldots, Y_n]'$ is the $n-1 \times 1 $ vector with the $i$-th row (of $\boldsymbol{Y}$) deleted; and $\mathbf{X}^{(-i)} = [\boldsymbol{X}_1, \ldots, \boldsymbol{X}_{i-1}, \boldsymbol{X}_{i+1}, \ldots, \boldsymbol{X}_n]'$ is the $(n-1)\times k$ matrix with the $i$-th row (of $\mathbf{X}$) deleted. Note that this is a very large system.
Once $\widehat{\boldsymbol{\beta}}^-$ is estimated (by least squares), the vector of residuals you want to estimate ($Y_i - \boldsymbol{X}_i'\widehat{\boldsymbol{\beta}}^{(-i)}$) can be written simply as $\mathrm{diag}(\mathbf{X}\mathrm{vec}(\widehat{\boldsymbol{\beta}}^-))$.
R code
Here is some R code to show how this can be done. The only trick here is constructing the augmented, row-deleted matrices, and then you are left with one (very) large least squares problem to solve.
iN <- 50 # number of observations
iK <- 4 # number of regressors (including constant)
mX <- matrix(rnorm(iN*iK), nrow = iN, ncol = iK) # design matrix
vBeta <- c(1, 2, 3, 4) # coefficients
vY <- mX%*%vBeta + matrix(rnorm(iN)) # dependent variable
mXAugmented <- (diag(iN)%x%mX)[-seq(from = 1,
to = iN*iN, by = iN+1), ] # augmented design matrix
vYAugmented <- vec(vY%*%t(rep(1, iN)))[-seq(from = 1,
to = iN*iN, by = iN+1), ] # augmented outcomes
vBetaMinus <- solve(t(mXAugmented)%*% mXAugmented,
t(mXAugmented)%*%vYAugmented) # estimated coefficients
mBetaMinus <- matrix(vBetaMinus, nrow = iK, ncol = iN)
vEpsilonAugmented <- diag(mX%*%mBetaMinus) # required residuals
matplot(vEpsilonAugmented, type = "l") # plot the estimated residuals