Asymptotic dist of an average involving OLS coefs? Suppose that we have iid sample of size $n$. i.e., the random vector $(Y_{i}, X_{1i}, X_{2i}, X_{3i})$ is iid from $1,\ldots,n$. And suppose the following relationship is true:
$$
Y_i = \beta_0 + \beta_1*X_{1i} + \beta_2*X_{2i} + \beta_3*X_{1i}*X_{2i} + \epsilon_i
$$
Suppose for simplicity that $X_{1i}$ and $X_{2i}$ are uniformly distributed from 0 to 1, and are correlated. Let's assume further that $\epsilon_i$ is normally distributed and independent of $X_{1i}$ and $X_{2i}$.
Let the OLS estimators be $\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2$.
Let $Z_i$ be
$$
Z_i = 1*\hat{\beta}_0 + 2*X_{1i}*\hat{\beta}_1 + 3*X_{2i}*\hat{\beta}_2 + 4*\hat{\beta}_3*X_{1i}*X_{2i}
$$
How do I find the asymptotic distribution of $\bar{Z}=\frac{1}{n}\sum_{i=1}^n Z_i$?
I cannot apply a CLT since the $Z_i$ are correlated with each other because of the $\hat{\beta}$. In addition to solving this particular case, any reference to theory I can study related to this would be helpful. I do not have an advanced statistical theory knowledge.
I would like to derive a non-degenerate asymptotic distribution, i.e., something like $\sqrt{n}(\bar{Z} - E(Z_i))$.
 A: The general form of the distribution of an individual quantity of this kind is quite complicated, and does not have a simple closed form.  I will give you an account of the general derivation and then show how the general form of the density.  Since the density is complicated, it is best to estimate it via simulation to compute the mean and variance, and then apply the central limit theorem to obtain an approximating distribution for the mean of the quantities.  I will generalise your problem in the initial analysis by allowing the explanatory variables to have any distribution.

General form of the distribution for an individual quantity: The simplest thing to do here is to start by conditioning on the explanatory variables and use standard regression results to find the conditional distribution of the quantity of interest.  You can then apply the law of total probability to find the marginal distribution.  To do this, note that when you formulate the OLS estimator you get the associated covariance matrix:
$$\mathbb{V}(\hat{\boldsymbol{\beta}}) = \sigma^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1}.$$
Let $\mathbf{x}_i^* \equiv \begin{bmatrix} 1 & 2 x_{1i} & 3 x_{2i} & 4 x_{1i} x_{2i} \end{bmatrix}$ denote the explanatory vector for the problem you are looking at.  You can write the random variable of interest as $Z_i = \mathbf{x}_i^* \hat{\boldsymbol{\beta}}$, which has the conditional distribution:
$$Z_i|\mathbf{x}_i^* \sim \mathcal{N} \Big( \mathbf{x}_i^* \boldsymbol{\beta}, \sigma^2 \mathbf{x}_i^* (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}_i^{* \text{T}} \Big).$$
The scalar form of these conditional moments is complicated, so we will leave things in matrix form.  Applying the law of total probability gives the marginal density function:
$$\begin{align}
f_{Z_i}(z) 
&= \int \mathcal{N} \Big( z \Big| \mathbf{x}_i^* \boldsymbol{\beta}, \sigma^2 \mathbf{x}_i^* (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}_i^{* \text{T}} \Big) dF (\mathbf{x}_i^*) \\[6pt]
&= \frac{1}{4 \pi^2} \int \det (\sigma^2 \mathbf{x}_i^* (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}_i^{* \text{T}} )^{-1/2} \exp \Bigg( - \frac{1}{2} \cdot \frac{(z - \mathbf{x}_i^* \boldsymbol{\beta})^2}{\sigma^2 \mathbf{x}_i^* (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}_i^{* \text{T}}} \Bigg) dF (\mathbf{x}_i^*). \\[6pt]
\end{align}$$
This gives a general equation for the density function that is a function of the underlying true regression parameters and the distribution of the explanatory variables.  The density does not have a closed form, including in the case you specify where the explanatory variables are standard uniform random variables.  The density function can be computed numerically from this equation, but it is simpler to proceed by simulation.  This can be done using the above conditional distribution, or it can be done by computing the intermediate OLS estimator and then applying the deterministic equation for the quantity of interest.

General form of the distribution for the mean of these quantities: Since your goal is to find the distribution of the standardised sample mean of these quantities, you can apply the central limit theorem.  Use simulations to estimate the true mean and variance from the underlying regression parameters and then apply the normal approximation to obtain the asymptotic distribution.  Given values for the vector xstar, the design matrix x, and the parameters beta and sigma you can simulate n values of $Z_i$ using the following function:
SIM <- function(n, xstar, x, beta, sigma) {
  MEAN <- sum(xstar*beta);
  VAR  <- sigma^2*(xstar %*% solve((t(x) %*% x), t(xstar)));
  rnorm(n, mean = MEAN, sd = sqrt(VAR)); }

