Analytically linking coefficients from alternative linear models (OLS) The general problem:
I have two alternative models I could use for my estimation
Model A: $y = \alpha^A+ X \beta^A_0 + Z\beta^A_1 + \varepsilon^A$
Model B: $y = \alpha^B + X \beta^B_0 + \varepsilon^B$
$X$: $(n\times k_X)$ and
${Z}$: $(n\times k_{Z})$
I want to estimate the effect of X on y, but can optionally use a set of  variables ($Z$) as control variables.
I am trying to assess what difference it makes for my slope estimate whether I estimate it through model A or model B. In particular I want to find an analytical expression to evaluate:


*

*The covariance between $\hat{\beta_0^A}$ and $\hat{\beta_0^B}$ 

*The conditional expectation: E[$\hat{\beta_0^A}$|$\hat{\beta_0^B}$]

*Or, ideally the full joint distributions / functional relationship between the two estimators


A simplified problem plus solution:
I already found a solution to the simplest case with $k_X=1$, and  $k_{Z}=1$, by moving from matrix notation to sum-notation:
OLS-slope coefficient estimate from model B (no control variables):
$$\beta^{B}_0 = \frac{\sum_ix_iy_i}{\sum_ix_i^2}$$
OLS-slope coefficient estimate from model A (with one control):
$$\beta^A_0 = \frac{(\sum_ix_iy_i)(\sum_iz_i^2)-(\sum_iy_iz_i)(\sum_iz_ix_i)}{(\sum_ix_i^2)(\sum_iz_i^2)-(\sum_ix_iz_i)^2}$$
where lower case letters indicate deviations from the mean, e.g. $x_i = X_i-\bar{X}$.
Rearranging terms yields:
$$\beta^A_0\left(1-\frac{(\sum_ix_iz_i)^2}{\sum_ix_i^2\sum_iz_i^2}\right) +  \frac{\sum_iy_iz_i\sum_iz_ix_i}{\sum_ix_i^2\sum_iz_i^2}= \beta^B_0
$$
This is what I want to generalize to the general case.
Edit: in the original version of the question I asked for the case where the two models have two alternative sets of covariates.
 A: To make things a little simpler notation-wise, let $X$ also contain a vector of ones (for the intercept). Also, let $\theta_1=(\alpha, \beta_0')'$ and $\theta=(\theta_1', \beta_1)'$ and $P_A=A(A'A)^{-1}A'$.
\begin{align*}
\hat\theta^A&=\begin{pmatrix}X'X & X'Z \\ Z'X & Z'Z\end{pmatrix}^{-1}\begin{pmatrix}X'y\\Z'y\end{pmatrix}\\
&=\begin{pmatrix}[X'X-X'Z(Z'Z)^{-1}Z'X]^{-1} & -[X'X-X'Z(Z'Z)^{-1}Z'X]^{-1}X'Z(Z'Z)^{-1}\\-[Z'Z-Z'X(X'X)^{-1}X'Z]^{-1}Z'X(X'X)^{-1} & [Z'Z-Z'X(X'X)^{-1}X'Z]^{-1}\end{pmatrix}\begin{pmatrix}X'y\\Z'y\end{pmatrix}\\
\end{align*}
So
\begin{align*}
\hat\theta^A_1&=[X'X-X'Z(Z'Z)^{-1}Z'X]^{-1}X'y -[X'X-X'Z(Z'Z)^{-1}Z'X]^{-1}X'Z(Z'Z)^{-1}Z'y\\
&=[X'X-X'Z(Z'Z)^{-1}Z'X]^{-1}(X'X)(X'X)^{-1}X'y -[X'X-X'Z(Z'Z)^{-1}Z'X]^{-1}X'Z(Z'Z)^{-1}Z'y\\
&=[X'(I-P_Z)X]^{-1}X'(I-P_Z)y\\
&=[X'(I-P_Z)X]^{-1}X'\left[X\hat\theta_1^B-P_Zy\right],
\end{align*}
where
\begin{align*}
\hat\theta_1^B=(X'X)^{-1}X'y.
\end{align*}
The covariance is then fairly simple:
$$
Cov(\hat\theta_1^A, \hat\theta_1^B)=Cov([X'(I-P_Z)X]^{-1}X'(I-P_Z)y, (X'X)^{-1}X'y)=[X'(I-P_Z)X]^{-1}X'(I-P_Z)V(y)X(X'X)^{-1}.
$$
Under the classical assumptions, $V(y)=\sigma^2I$, so
\begin{align*}
Cov(\hat\theta_1^A, \hat\theta_1^B)&=\sigma^2[X'(I-P_Z)X]^{-1}X'(I-P_Z)X(X'X)^{-1}=\sigma^2(X'X)^{-1}.
\end{align*}
When $X$ and $Z$ are both vectors, we get
\begin{align*}
\hat\beta_0^A&=\left(\sum_{i=1}^nx_i^2-\frac{(\sum_{i=1}^nx_iz_i)^2}{\sum_{i=1}^nz_i^2}\right)^{-1}\left(\hat\beta_0^B\sum_{i=1}x_i^2-\frac{\sum_{i=1}^nx_iz_i\sum_{i=1}^nz_iy_i}{\sum_{i=1}^nz_i^2}\right)\\
\hat\beta_0^A&=\left(1-\frac{(\sum_{i=1}^nx_iz_i)^2}{\sum_{i=1}^nx_i^2\sum_{i=1}^nz_i^2}\right)^{-1}\left(\hat\beta_0^B-\frac{\sum_{i=1}^nx_iz_i\sum_{i=1}^nz_iy_i}{\sum_{i=1}^nx_i^2\sum_{i=1}^nz_i^2}\right)\\
\hat\beta_0^A\left(1-\frac{(\sum_{i=1}^nx_iz_i)^2}{\sum_{i=1}^nx_i^2\sum_{i=1}^nz_i^2}\right)&=\hat\beta_0^B-\frac{\sum_{i=1}^nx_iz_i\sum_{i=1}^nz_iy_i}{\sum_{i=1}^nx_i^2\sum_{i=1}^nz_i^2}\\
\hat\beta_0^A\left(1-\frac{(\sum_{i=1}^nx_iz_i)^2}{\sum_{i=1}^nx_i^2\sum_{i=1}^nz_i^2}\right)+\frac{\sum_{i=1}^nx_iz_i\sum_{i=1}^nz_iy_i}{\sum_{i=1}^nx_i^2\sum_{i=1}^nz_i^2}&=\hat\beta_0^B,
\end{align*}
which is what you got for the same case.
