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.
