Although your question asks about a model where the error terms are independent, I'm going to generalise to drop this assumption in order to give a slightly broader answer. This will allow me to show you the effect of that assumption on your analysis more easily, and I hope it will also give a bit more insight into the interactions occurring in the multivariate regression model.
Also, at the outset it is worth noting that the analytical results for the multivariate regression model turn out to be mostly identical to —or analogous to— the case where we have a single response variable. This is true of the OLS estimator, the hat matrix, the error variance estimator, and many of the other statistical results. I'll show some of this below.
For your problem, let's start by writing your overall model equation in matrix form:
$$\begin{bmatrix} y_{1,1} & y_{1,2} \\ y_{2,1} & y_{2,2} \\ \vdots & \vdots \\ y_{n,1} & y_{n,2} \end{bmatrix}
= \begin{bmatrix} 1 & x_1 \\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \end{bmatrix} \begin{bmatrix} a & c \\ b & d \end{bmatrix}
+ \begin{bmatrix} \epsilon_{1,1} & \epsilon_{1,2} \\ \epsilon_{2,1} & \epsilon_{2,2} \\ \vdots & \vdots \\ \epsilon_{n,1} & \epsilon_{n,2} \end{bmatrix}
\quad \quad \quad \quad \quad
\begin{bmatrix} \epsilon_{i,1} \\ \epsilon_{i,2} \end{bmatrix} \sim \text{IID N}(\mathbf{0}, \boldsymbol{\Sigma}).$$
Here I am using a general variance matrix $\boldsymbol{\Sigma}$ for the error variance, without assuming that the errors of the two parts are uncorrelated (but we can look at the specific case of interest to you later). Defining the design matrix $\mathbf{x}$ and appropriate vectors for the responses, coefficients and errors of the two parts, we can write this model equation in more succinct form as:
$$\begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 \end{bmatrix}
= \mathbf{x} \begin{bmatrix} \boldsymbol{\beta}_1 & \boldsymbol{\beta}_2 \end{bmatrix}
+ \begin{bmatrix} \boldsymbol{\epsilon}_1 & \boldsymbol{\epsilon}_2 \end{bmatrix}.$$
In multivariate regression, it turns out that the OLS estimator for the overall model is the same as imposing OLS estimation on the two parts. The OLS estimator is given by:$^\dagger$
$$\begin{align}
\begin{bmatrix} \hat{\boldsymbol{\beta}}_1 & \hat{\boldsymbol{\beta}}_2 \end{bmatrix}
&= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T}
\begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 \end{bmatrix} \\[6pt]
&= \begin{bmatrix} (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \mathbf{y}_1 & (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T} \mathbf{y}_2 \end{bmatrix}. \\[6pt]
\end{align}$$
Assuming the model is correctly specified, the OLS estimator is an unbiased estimator of the true coefficients. If we write the coefficient estimators as a single vector $\hat{\boldsymbol{\beta}} = [\hat{a},\hat{b},\hat{c},\hat{d}]^\text{T}$, this vector has variance matrix:
$$\begin{align}
\mathbb{V}(\hat{\boldsymbol{\beta}})
&= (\mathbf{x}^\text{T} \mathbf{x})^{-1} \otimes \boldsymbol{\Sigma}
= \begin{bmatrix} \sigma_1^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} & \rho \sigma_1 \sigma_2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \\ \rho \sigma_1 \sigma_2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} & \sigma_2^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \end{bmatrix}.
\end{align}$$
Moreover, under the Gaussian version of the model the OLS estimator is normally distributed, so we have:
$$\hat{\boldsymbol{\beta}} - \boldsymbol{\beta}
\sim \text{N} \Bigg( \begin{bmatrix} \mathbf{0} \\[4pt] \mathbf{0} \end{bmatrix}, \begin{bmatrix} \sigma_1^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} & \rho \sigma_1 \sigma_2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \\ \rho \sigma_1 \sigma_2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} & \sigma_2^2 (\mathbf{x}^\text{T} \mathbf{x})^{-1} \end{bmatrix} \Bigg).$$
From the above form, we can see that the estimators for the two coefficient vectors in the model parts have covariance $\mathbb{C}(\hat{\boldsymbol{\beta}}_1, \hat{\boldsymbol{\beta}}_2) = \rho \sigma_{1} \sigma_{2} (\mathbf{x}^\text{T} \mathbf{x})^{-1}$. Consequently, in the case where the underlying errors for the two model parts are uncorrelated (i.e., when $\rho = 0$) the coefficient estimators for the two model parts are also uncorrelated.
In order to form a confidence set for $\boldsymbol{\beta}$ you also need to estimate the error variance matrix $\boldsymbol{\Sigma}$. The standard estimator is the adjusted MLE (adjusted to be unbiased), which is:
$$\hat{\boldsymbol{\Sigma}}
= \frac{\begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 \end{bmatrix}^\text{T} (\mathbf{I}_n - \mathbf{x} (\mathbf{x}^\text{T} \mathbf{x})^{-1} \mathbf{x}^\text{T}) \begin{bmatrix} \mathbf{y}_1 & \mathbf{y}_2 \end{bmatrix}}{n-2}.$$
Using these results you can form a confidence set for $\boldsymbol{\beta}$, which will be a hyper-ellipse in four-dimensional Euclidean space.
$^\dagger$ It's also worth noting that the OLS estimator is equivalent to the MLE in the Gaussian case, just as in standard regression.