# Bootstrapping multivariate multiple variable regression

I would like to know how to bootstrap multiple variable and multivariate regression. Let me first go over the single variable case. I realize that in general either the (x,y) pair can be bootstrapped, or only the error term. I will focus on the case for bootstrapping the error term, as that is in line with Gauss Markov conditions.

$y = \beta x_1 + \epsilon$

Where y is the dependent variable, x is the independent variable, epsilon is the error term and beta is the coefficient, which we estimate. After estimating beta, we can re-sample the residuals with replacement. Let's denote the re-sampled residual as $eB$. Now we can estimate a new equation:

$(y - eB) = \beta x_1 + \epsilon$

The equation most likely will yield different beta coefficient than the original one. We can repeat the process of re-sampling and estimating the equation multiple (thousands) of times to derive properties such as confidence interval. Very effective if the data is not believed to be normally distributed.

Now, for the multiple variable case. Which is characterized, for example by the following equation:

$y = \beta_0+ \beta_1 x_1 + \beta_2 x_{2}+ \epsilon$

Again, we can re-sample the error terms without any issues. After re-sapmling, the re-sampled error terms can again be inserted into the equation. We have:

$(y-eB) = \beta_0+ \beta_1 x_1 + \beta_2 x_{2}+ \epsilon$

Now, if we want to find the confidence interval for $\beta_1$. how do we proceed? Is it the case that we simply hold all the other coefficients fixed and re-estimate only $\beta_1$ using the re-sampled error terms? Use the re-sampled errors again and again, to do the same for other coefficients and for the coefficients as a whole. Rinse and repeat till you have enough estimates to form a confidence interval. Is this correct?

Now, for the multivariate case. As an illustration let's use the following VAR equation:

$$\begin{cases} &y_{1t} = \alpha_0+ \alpha_1 y_{1,t-1} + \alpha_2 y_{2,t-1}+ \epsilon_{1t} &\\\ &y_{2t} = \alpha_5+ \alpha_3 y_{1,t-1} + \alpha_4 y_{2,t-1}+ \epsilon_{2t} \end{cases}$$

Is it the case that in this situation both of the error terms have to be re-sampled. Afterwards, proceed as in the multiple variable case for both equations independently? This, even if we orthogonalize the error terms (using the Cholesky decomposition, for example).