# Linear model with random coefficient

Suppose I wanted to fit a model of the form $$y_i = \beta_0 + \sum_{1 \le j \le k} \beta_j X_{i,j} + \gamma_i Z_i + \epsilon_i,$$ to some data, where the regressors $X$ and $Z$, and the regressand $y$ are observed, and where $\gamma_i$ is a Bernoulli random variable that equals one with (unknown) probability $p$ and is zero otherwise. We can assume all kinds of 'regularity': the errors $\epsilon_i$ are independent of the regressors, and are independent of $\gamma_i$, etc.

Some questions:

1. What is the name for this, if there is one?
2. If I just throw the $Z$ data in with the $X$ data and perform an ordinary multiple least squares, will the least squares coefficient corresponding to the $Z$ term converge to $p$ as I add more observations?
3. If least squares regression is advisable for this model, what is the distribution of the least squares coefficient corresponding to the $Z$ term under the null hypothesis $p = p_0$? (For a 'deterministic' regression, the coefficient has a certain t-distribution, with parameters depending on $p_0, n, k$ and the design matrix; I am looking for the analogue for the random coefficient.)
4. If least squares regression is advisable, how will the presence of the random coefficient affect the distribution of the other (sample) regression coefficients?
-

I will try to answer the questions 2 to 4. Suppose that we observe sample $(y_i,\mathbf{x}_i,z_i,\gamma_i,\varepsilon_i)$. Suppose that our model is

$$y_i=\mathbf{x}_i\beta+\gamma_iz_i+\varepsilon_i$$

and

$$E(\varepsilon_i|\mathbf{x}_i,z_i,\gamma_i)=0.$$

The least squares estimate of the regression will be \begin{align} (\hat{\beta},\hat{\gamma})'=\left(\sum_{i=1}^n \begin{bmatrix} \mathbf{x}_i'\\ z_i \end{bmatrix} [\mathbf{x}_i,z_i]\right)^{-1}\sum_{i=1}^n\begin{bmatrix} \mathbf{x}_i'\\ z_i \end{bmatrix}y_i \end{align}

Now since we assumed sample due to law of large numbers we get that

\begin{align} \frac{1}{n}\sum_{i=1}^n \begin{bmatrix} \mathbf{x}_i'\\ z_i \end{bmatrix} [\mathbf{x}_i,z_i]\to \begin{bmatrix} E\mathbf{x}_1'\mathbf{x}_1 & E\mathbf{x_1}'z_1\\ E\mathbf{x}_1z_1 & Ez_1^2 \end{bmatrix} \end{align}

Now \begin{align} \sum_{i=1}^n\begin{bmatrix} \mathbf{x}_i'\\ z_i \end{bmatrix}y_i=\sum_{i=1}^n\begin{bmatrix} \mathbf{x}_i'\\ z_i \end{bmatrix}\varepsilon_i+\sum_{i=1}^n \begin{bmatrix} \mathbf{x}_i'\mathbf{x}_i\beta +\mathbf{x}_i'z_i\gamma_i \\ \mathbf{x}_iz_i\beta + z_i^2\gamma_i \end{bmatrix} \end{align}

Due to law of large numbers and our conditional expectation condition we get that \begin{align} \frac{1}{n}\sum_{i=1}^n\begin{bmatrix} \mathbf{x}_i'\\ z_i \end{bmatrix}\varepsilon_i\to 0. \end{align}

Now comes the part where we need more assumptions. Assume that $(\mathbf{x}_i,z_i)$ is independent of the $\gamma_i$. Then due to law of large numbers

\begin{align} \frac{1}{n}\sum_{i=1}^n\mathbf{x}_iz_i\gamma_i\to E\mathbf{x}_1z_1\gamma_1=pE\mathbf{x}_1z_1 \end{align}

where $p=P(\gamma_i=1)$. Similarly

\begin{align} \frac{1}{n}\sum_{i=1}^nz_i^2\gamma_i\to Ez_1^2\gamma_1=pEz_1^2 \end{align}

Gathering all the results we get

\begin{align} (\hat{\beta},\hat{\gamma})'\to \left(\begin{bmatrix} E\mathbf{x}_1'\mathbf{x}_1 & E\mathbf{x_1}'z_1\\ E\mathbf{x}_1z_1 & Ez_1^2 \end{bmatrix}\right)^{-1}\begin{bmatrix} E\mathbf{x}_1'\mathbf{x}_1\beta + pE\mathbf{x_1}'z_1\\ E\mathbf{x}_1z_1\beta + pEz_1^2 \end{bmatrix}=(\beta,p)' \end{align}

So the answer to second question is yes. Simple experiment in R confirms this:

> g<-sample(0:1,1000,prob=c(1/3,1-1/3),replace=TRUE)
> z <- rnorm(1000)
> y<-1+g*z+rnorm(1000)/3
> dt <- data.frame(y=y,z=z)
> lm(y~z,data=dt)

Call:
lm(formula = y ~ z, data = dt)

Coefficients:
(Intercept)            z
1.0050       0.6595


Now let us proceed to question 3. Introduce notation $X_i=(\mathbf{x}_i,z_i)$. Then

\begin{align} \sqrt{n}\begin{bmatrix} \hat{\beta}-\beta\\ \hat{\gamma}-p \end{bmatrix}=\left(\frac{1}{n}\sum_{i=1}^nX_i'X_i\right)^{-1}\frac{1}{\sqrt{n}}\left(\sum_{i=1}^nX_i'\varepsilon_i+\sum_{i=1}^n \begin{bmatrix} \mathbf{x}_i'z_i\\\ z_i^2 \end{bmatrix}(\gamma_i-p)\right) \end{align}

Introduce notation $Z_i=(\mathbf{x}_iz_i,z_i^2)$ and $C_i=(X_i\varepsilon_i,Z_i(\gamma_i-p))$. We have that $EC_i=0$ and since we have sample we can apply multivariate central limit theorem for $C_i$:

\begin{align} \frac{1}{\sqrt{n}}\sum_{i=1}^nC_i'\to N(0,\Sigma_C) \end{align} where \begin{align} \Sigma_C=EC_i'C_i&= \begin{bmatrix} EX_1'X_1\varepsilon_1^2 & EX_1'Z_1\varepsilon_1(\gamma_1-p)\\ EZ_1'X_1\varepsilon_1(\gamma_1-p) & EZ_1'Z_1(\gamma_1-p)^2 \end{bmatrix}\\ &=\begin{bmatrix} EX_1'X_1\varepsilon_1^2 & 0\\ 0 & EZ_1'Z_1(\gamma_1-p)^2 \end{bmatrix} \end{align} due to condition $E(\varepsilon_i|\mathbf{x}_i,z_i,\gamma_i)=0$.

Now assume further that $E(\varepsilon_i^2|\mathbf{x}_i,z_i,\gamma_i)=\sigma^2$. Denote \begin{align} \mathbf{A}=EX_1'X_1, \quad \mathbf{B}=EZ_1'Z_1 \end{align} Then we get \begin{align} \left(\frac{1}{n}\sum_{i=1}^nX_i'X_i\right)^{-1}\frac{1}{\sqrt{n}}\left(\sum_{i=1}^nX_i'\varepsilon_i+\sum_{i=1}^n \begin{bmatrix} \mathbf{x}_i'z_i\\\ z_i^2 \end{bmatrix}(\gamma_i-p)\right)\to N(0,\sigma^2A+p(1-p)B) \end{align} so \begin{align} \sqrt{n}\begin{bmatrix} \hat{\beta}-\beta\\ \hat{\gamma}-p \end{bmatrix}\to N(0,\sigma^2A^{-1}+p(1-p)A^{-1}BA^{-1}) \end{align}

This should give an idea how to construct the test for testing $H_0:p=p_0$. I am not happy about the presence of $p$ in the variance matrix, this could pose some problems.

Now given the above the answer to the fourth question is that random coefficient influences only the covariance matrix of the other coefficients. If there is bias, it vanishes asymptotically.

I should note that the way I derived these results is pretty straightforward application of LLN and CLT. If there is some elegant way to avoid this I would really like to know.

-