Why can't I use OLS for switching regression? I am not sure about the exact reason, why I can't use the simple OLS for switching regression.
I have two outcome equations:
$Y_1=\beta_0 + \beta_1'X + U_1$
$Y_0=\alpha_0 + \alpha_1'X+U_0$
and a selection equation:
$D=1$ if $\gamma'Z+\epsilon>0$ and zero otherwise.
Why can't I use simple OLS? The correct approach is to use a 2 step approach or a maximum likelihood (ML) approach.
Is it, because I have a selection bias and therefore the $U$ are correlated with the $X$?
 A: Simply put: the simple OLS estimator will be inconsistent when $U_{1}$
or $U_{0}$ is not mean indepent with $\varepsilon$, i..e $\mathrm{E}\left(U_{1}\mid\varepsilon\right)\neq0$
or $\mathrm{E}\left(U_{0}\mid\varepsilon\right)\neq0$.
To give more detailed explanation, let me change your notations slightly
for expositional simplicity. Denote the two outcome equations:
$$
Y_{1}=X'\beta+U_{1},
$$
$$
Y_{0}=X'\alpha+U_{0}.
$$
Also assume $X$ is independent of $U_{1}$, $U_{0}$, $\varepsilon$
and $Z$. Denote $i$ a unit of observation. Using the selection variable
$D$, we have one single regression
$$
Y_{i}=D_{i}X_{i}'\beta+\left(1-D_{i}\right)X_{i}'\alpha+D_{i}U_{1}+\left(1-D_{i}\right)U_{0}.
$$
The simple OLS minimizes
$$
\sum_{i=1}^{N}\left[Y_{i}-\left\{ D_{i}X_{I}'\beta+\left(1-D_{i}\right)X_{l}'\alpha\right\} \right]^{2}.
$$
The minimizer or the simple OLS estimator of $\lambda'=\left(\beta,\alpha\right)'$
is
$$
\lambda_{N}=\left(\mathbf{W}'\mathbf{W}\right)^{-1}\mathbf{W}'\mathbf{Y},
$$
where $\mathbf{Y}=\left(Y_{1},\ldots,Y_{N}\right)'$ and $\mathbf{W}$
is a $N\times(2K)$ matrix, $K$ is the dimension of $\alpha$ or
$\beta$,
$$
\mathbf{W}=\left[\begin{array}{cc}
D_{1}X_{1}' & \left(1-D_{1}\right)X_{1}'\\
\vdots & \vdots\\
D_{N}X_{N}' & \left(1-D_{N}\right)X_{N}'
\end{array}\right].
$$
Denote $V_{i}=D_{i}U_{1,i}+\left(1-D_{i}\right)U_{0,i}$, and $\mathbf{V}'=\left(V_{1},\ldots,V_{N}\right)'$.
The estimator $\lambda_{N}$ can be written as follows,
$$
\left(\mathbf{W}'\mathbf{W}\right)^{-1}\mathbf{W}'\mathbf{Y}=\lambda + \left(\frac{1}{N}\mathbf{W}'\mathbf{W}\right)^{-1}\left(\frac{1}{N}\mathbf{W}'\mathbf{V}\right),
$$
since $\mathbf{Y}=\mathbf{W}'\lambda+\mathbf{V}$. Denote $\mathrm{plim}_{N\rightarrow\infty}\mathbf{W}'\mathbf{W}/N=\mathbf{A}$.
By the law of large numbers,
$$
\lambda_{N}\rightarrow_{p}\mathbf{A}^{-1}\mathrm{{E}}\left(\left[\begin{array}{c}
D_{i}X_{i}\left(D_{i}U_{1,i}+\left(1-D_{i}\right)U_{0,i}\right)\\
\left(1-D_{i}\right)X_{i}\left(D_{i}U_{1,i}+\left(1-D_{i}\right)U_{0,i}\right)
\end{array}\right]\right).
$$
By the law of itereated expectation, we have
\begin{eqnarray*}
\mathrm{{E}}\left\{ D_{i}X_{i}\left(D_{i}U_{1,i}+\left(1-D_{i}\right)U_{0,i}\right)\right\}  & = & \mathrm{E}_{D_{i}}\left[\mathrm{{E}}\left\{ D_{i}X_{i}\left(D_{i}U_{1,i}+\left(1-D_{i}\right)U_{0,i}\right)\mid D_{i} \right\} \right]\\
 & = & \Pr\left(D_{i}=1\right)\mathrm{{E}}\left(X_{i}U_{1,i}\mid D_{i}=1\right)\\
 & = & \Pr\left(D_{i}=1\right)\mathrm{{E}}\left(X_{i}\right)\mathrm{{E}}\left(U_{1,i}\mid D_{i}=1\right).
\end{eqnarray*}
The inconsistency comes from the fact that
$$
\mathrm{{E}}\left(U_{1,i}\mid D_{i}=1\right)=\mathrm{{E}}\left(U_{1,i}\mid\varepsilon>-\gamma'Z_{i}\right)\neq0
$$
in general when $U_{1}$ and $\varepsilon$ are correlated. For example,
suppose $\left(U_{1},U_{0},\varepsilon\right)'$ follows a trivariate
normal distribution with means $0$ and covariance matrix being
$$
\left[\begin{array}{ccc}
\sigma_{1}^{2} & \sigma_{12} & \sigma_{1\varepsilon}\\
 & \sigma_{0}^{2} & \sigma_{0\varepsilon}\\
 &  & 1
\end{array}\right].
$$
We have $\mathrm{{E}}\left(U_{1,i}\mid\varepsilon>-\gamma'Z_{i}\right)=-\sigma_{1\varepsilon}\phi\left(-\gamma'Z_{i}\right)/\Phi\left(-\gamma'Z_{i}\right)$.
Here $\phi$ and $\Phi$ are the pdf and CDF of a standard normal.
Similar arguments can be applied to show $\mathrm{{E}}\left\{ \left(1-D_{i}\right)X_{i}\left(D_{i}U_{1,i}-\left(1-D_{i}\right)U_{0,i}\right)\right\} \neq0$.
