[proofplan]
We use the column-stacked representation of the SUR system, where the common-regressor design is $Z=I_p\otimes X$ and the GLS weight is $\Sigma^{-1}\otimes I_n$. The key computation is that the GLS normal matrix factors as $\Sigma^{-1}\otimes X^\top X$, while the GLS right-hand side factors as $(\Sigma^{-1}\otimes X^\top)\operatorname{vec}(Y)$. Inverting the Kronecker product cancels the factor $\Sigma^{-1}$, leaving precisely the vectorization of the equation-by-equation OLS matrix $(X^\top X)^{-1}X^\top Y$.
[/proofplan]
[step:Verify that all matrices in the GLS formula are well defined]
Since $X$ has full column rank, the Gram matrix $X^\top X \in \mathbb{R}^{k \times k}$ is invertible. Since $\Sigma$ is symmetric positive definite, $\Sigma^{-1} \in \mathbb{R}^{p \times p}$ exists.
The stacked design matrix $Z=I_p\otimes X$ has size $np \times kp$, and the GLS weight matrix $\Sigma^{-1}\otimes I_n$ has size $np \times np$. Therefore
\begin{align*}
Z^\top(\Sigma^{-1}\otimes I_n)Z \in \mathbb{R}^{kp \times kp}
\end{align*}
and
\begin{align*}
Z^\top(\Sigma^{-1}\otimes I_n)\operatorname{vec}(Y) \in \mathbb{R}^{kp}.
\end{align*}
We will compute the normal matrix explicitly and see that it is invertible.
[/step]
[step:Factor the GLS normal matrix into equation and regressor components]
Using $Z=I_p\otimes X$ and $(I_p\otimes X)^\top=I_p\otimes X^\top$, we compute
\begin{align*}
Z^\top(\Sigma^{-1}\otimes I_n)Z
&=(I_p\otimes X^\top)(\Sigma^{-1}\otimes I_n)(I_p\otimes X) \\
&=\Sigma^{-1}\otimes X^\top X.
\end{align*}
The second equality follows from the defining multiplication rule for Kronecker products:
\begin{align*}
(A\otimes B)(C\otimes D)=(AC)\otimes(BD),
\end{align*}
applied first with $A=I_p$, $B=X^\top$, $C=\Sigma^{-1}$, $D=I_n$, and then with $A=\Sigma^{-1}$, $B=X^\top$, $C=I_p$, $D=X$.
Because both $\Sigma^{-1}$ and $X^\top X$ are invertible, the Kronecker product $\Sigma^{-1}\otimes X^\top X$ is invertible, with inverse
\begin{align*}
(\Sigma^{-1}\otimes X^\top X)^{-1}
=
\Sigma\otimes (X^\top X)^{-1}.
\end{align*}
[/step]
[step:Factor the GLS right-hand side]
Again using $Z=I_p\otimes X$, we obtain
\begin{align*}
Z^\top(\Sigma^{-1}\otimes I_n)\operatorname{vec}(Y)
&=(I_p\otimes X^\top)(\Sigma^{-1}\otimes I_n)\operatorname{vec}(Y) \\
&=(\Sigma^{-1}\otimes X^\top)\operatorname{vec}(Y).
\end{align*}
Thus the GLS coefficient vector is
\begin{align*}
\hat{\beta}_{\mathrm{GLS}}
&=
\left(\Sigma^{-1}\otimes X^\top X\right)^{-1}
(\Sigma^{-1}\otimes X^\top)\operatorname{vec}(Y) \\
&=
\left(\Sigma\otimes (X^\top X)^{-1}\right)
(\Sigma^{-1}\otimes X^\top)\operatorname{vec}(Y).
\end{align*}
[/step]
[step:Cancel the covariance weighting in the coefficient equation]
Apply the Kronecker product multiplication rule once more:
\begin{align*}
\left(\Sigma\otimes (X^\top X)^{-1}\right)
(\Sigma^{-1}\otimes X^\top)
&=
(\Sigma\Sigma^{-1})\otimes\left((X^\top X)^{-1}X^\top\right) \\
&=
I_p\otimes\left((X^\top X)^{-1}X^\top\right).
\end{align*}
Therefore
\begin{align*}
\hat{\beta}_{\mathrm{GLS}}
=
\left(I_p\otimes (X^\top X)^{-1}X^\top\right)\operatorname{vec}(Y).
\end{align*}
[/step]
[step:Identify the resulting vector as the vectorized OLS coefficient matrix]
Define the matrix $A \in \mathbb{R}^{k \times n}$ by
\begin{align*}
A := (X^\top X)^{-1}X^\top.
\end{align*}
For every matrix $Y \in \mathbb{R}^{n \times p}$, left multiplication by $A$ acts column by column under column-stacking:
\begin{align*}
\operatorname{vec}(AY)=(I_p\otimes A)\operatorname{vec}(Y).
\end{align*}
Applying this identity with $A=(X^\top X)^{-1}X^\top$ gives
\begin{align*}
\hat{\beta}_{\mathrm{GLS}}
&=
\left(I_p\otimes (X^\top X)^{-1}X^\top\right)\operatorname{vec}(Y) \\
&=
\operatorname{vec}\left((X^\top X)^{-1}X^\top Y\right).
\end{align*}
By definition, $\operatorname{vec}(\hat B_{\mathrm{GLS}})=\hat{\beta}_{\mathrm{GLS}}$. Since the column-stacking map is injective, the equality of vectorizations implies
\begin{align*}
\hat B_{\mathrm{GLS}}=(X^\top X)^{-1}X^\top Y.
\end{align*}
This is exactly the equation-by-equation ordinary least squares estimator for the common regressor matrix $X$.
[/step]
