The restriction just says that the variables in $X_1$ do not enter the regression, so that the formula constrains those coefficients to zero, and the others then are the standard OLS coefficients of a regression of $y$ on $X_2$.
Here is an implementation in R (a more formal answer follows below):
library(lrmest)
# example data
n <- 100
k_1 <- 3
k_2 <- 4
X_1 <- matrix(rnorm(n*k_1),ncol=k_1)
X_2 <- matrix(rnorm(n*k_2),ncol=k_2)
y <- rnorm(n)
X <- cbind(X_1,X_2)
# implementing the restriction in lm directly
coef(lm(y~X_2-1))
>
X_21 X_22 X_23 X_24
-0.02171519 -0.04595167 0.08971665 0.04871940
# restricted least squares
R <- cbind(diag(k_1),matrix(0,nrow=k_1,ncol=k_2))
r <- delta <- rep(0,k_1)
rls(y~X-1,r,R,delt=delta)
>
$`*****Restricted Least Square Estimator*****`
Estimate Standard_error t_statistic pvalue
X1 0.0000 0.0000 0 1
X2 0.0000 0.0000 0 1
X3 0.0000 0.0000 0 1
X4 -0.0217 0.0960 NA NA
X5 -0.0460 0.1091 NA NA
X6 0.0897 0.1093 NA NA
X7 0.0487 0.0922 NA NA
Both question and answer make use of properties of restricted OLS, which are for example discussed further here.
In particular, the restricted OLS estimator is related to the unrestricted one via
$$
\beta^*=\hat\beta+(X^TX)^{-1}R^T[R(X^TX)^{-1}R^T]^{-1}(r-R\hat\beta)
$$
Algebraically, the mistake seems to be in
$$R\hat\beta - r = (I_l:0)\hat\beta,$$ which is just $\hat\beta_l$.
In more detail, denote
\begin{align*}
(X^TX)^{-1}&=\left( \begin{array} {c,c} X_1^TX_1&X_1^TX_2 \\ X_2^TX_1&X_2^TX_2\end{array} \right)^{-1}\\&\equiv\left( \begin{array} {c,c} \tilde A&\tilde B \\ \tilde C&\tilde D\end{array} \right)
\end{align*}
Then,
$$R(X^TX)^{-1}R^T = \tilde A$$
and
$$(X^TX)^{-1}R^T = \left( \begin{array} {c} \tilde A \\ \tilde C\end{array} \right)$$
leading to
\begin{align*}
\beta^*&= \left( \begin{array} {c} \hat\beta_l \\ \hat\beta_k\end{array} \right) - \left( \begin{array} {c,c} \tilde A \\ \tilde C\end{array} \right)\tilde A^{-1}R\hat\beta\\
&= \left( \begin{array} {c} \hat\beta_l \\ \hat\beta_k\end{array} \right) - \left( \begin{array} {c,c} I \\ \tilde C\tilde A^{-1}\end{array} \right)\hat\beta_l\\
\end{align*}
Thus, using
$$
\hat\beta=\left( \begin{array} {c,c} \tilde A&\tilde B \\ \tilde C&\tilde D\end{array} \right)\left( \begin{array} {c} X_1^Ty\\X_2^Ty\end{array} \right),
$$
we have
$$
\hat\beta_k=\tilde C X_1^Ty+\tilde DX_2^Ty
$$
and
$$
\hat\beta_l=\tilde A X_1^Ty+\tilde BX_2^Ty
$$
so that
$$
\tilde C\tilde A^{-1}\hat\beta_l=\tilde C\tilde A^{-1}\tilde A X_1^Ty+\tilde C\tilde A^{-1}\tilde BX_2^Ty
$$
and
\begin{align*}
\hat\beta_k-\tilde C\tilde A^{-1}\hat\beta_l&=\tilde C X_1^Ty+\tilde DX_2^Ty-
\tilde C\tilde A^{-1}\tilde A X_1^Ty-\tilde C\tilde A^{-1}\tilde BX_2^Ty\\
&=(\tilde D-\tilde C\tilde A^{-1}\tilde B)X_2^Ty
\end{align*}
Now, use results for partitioned inverses
\begin{align*}
(X^TX)^{-1}
&=\left( \begin{array} {c,c} \tilde A&\tilde B \\ \tilde C&\tilde D\end{array} \right)\\&=
\left( \begin{array} {c,c} A&B \\ C&D\end{array} \right)^{-1}\\
&=\left( \begin{array} {c,c} E&-EBD^{-1} \\ -D^{-1}CE&D^{-1}+D^{-1}CEBD^{-1}\end{array} \right)\\
\end{align*}
where $E=(A-BD^{-1}C)^{-1}$.
Hence,
\begin{align*}
\tilde D-\tilde C\tilde A^{-1}\tilde B&=D^{-1}+D^{-1}CEBD^{-1}-D^{-1}CEE^{-1}EBD^{-1}\\&=D^{-1}
\end{align*}
The result is complete by noting that
$$D^{-1}=(X_2^TX_2)^{-1}$$
so that
$$\hat\beta_k-\tilde C\tilde A^{-1}\hat\beta_l=(X_2^TX_2)^{-1}X_2^Ty$$