Proof Verification: Joint variance of the product of a random matrix with a random vector BACKGROUND
QUESTIONS
Is the proof of my claim correct?
How might my proof be improved?

Claim:   (1) The joint-covariance matrix of the product of a real random matrix $X$ of dimension $v\times m$ and a real random matrix $Y$ of dimension $m\times 1$ is a real matrix of dimension $v\times v$.            (2) The element   on the $k^\textrm{th}$ row and $l^\textrm{th}$ column of the joint-covariance matrix, which I denote as $\operatorname {E} \left[(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]   )(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]  )^{\top }\right]_{k,l}$, is given as
  $$\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\Bigl(
\operatorname {cov}_X( X_{ki},    X_{lj})
 + \operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
 \Bigr)\Bigl(
 \operatorname {cov}_Y( Y_{i},    Y_{j}
 )
 + \operatorname {E}_Y \left[ Y_{i}     \right]
\operatorname {E}_Y \left[ Y_{j}    \right]
 \Bigr) -\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right] \operatorname {E}_Y \left[ Y_{i}     \right]
\operatorname {E}_Y \left[ Y_{j}    \right]
$$

PROOF
PART I
By $\mathbf{X}$ I denote a real random matrix of dimension $v\times m$. By $\mathbf{Y}$ I denote a real random matrix of dimension $m\times 1$.  I write these matrices explicitilty as
\begin{align*}
\mathbf{X} 
&=
\begin{bmatrix}
X_{11} & \cdots & X_{1m}
\\
\vdots & \vdots & \vdots
\\
X_{v1} & \cdots & X_{vm}
\end{bmatrix},~\textrm{and}
\\
\mathbf{Y}
&=
\begin{bmatrix}
Y_{1} 
\\
\vdots  
\\
Y_{m}  
\end{bmatrix}~\textrm{respectively.}
\end{align*}
Apriori, I state that $X_{ij}$ and $Y_k$ are statistically independent for any and all $i$ in $1,\ldots, v$;  any and all $j$ in $1,\ldots, m$; and any and all $k$ in $1,\ldots, m$.
The product $\mathbf{X} \,\mathbf{Y}$ can be written explicitly as
\begin{align*}
\mathbf{X} \,\mathbf{Y}
&=
\begin{bmatrix}
X_{11} & \cdots & X_{1m}
\\
\vdots & \vdots & \vdots
\\
X_{v1} & \cdots & X_{vm}
\end{bmatrix}
\begin{bmatrix}
Y_{1} 
\\
\vdots  
\\
Y_{m}  
\end{bmatrix}
\\
&=
\begin{bmatrix}
\sum\limits_{i=1}^m X_{1i}\,Y_{i}
\\
\vdots  
\\
\sum\limits_{i=1}^m X_{vi}\,Y_{i}
\end{bmatrix}
\end{align*}
Adapting from [1], since $\mathbf{X} \,\mathbf{Y}$ is a vector-valued random vector, with values in $\mathbb{R}^v$, then a natural generalization of variance is 
$$ \operatorname {E} \left[(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]   )(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]  )^{\top }\right].
$$
From [2], since   $X_{ij}$ and $Y_k$ are statistically independent, 
\begin{align*}
\operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right] 
=&
\operatorname {E}_X \left[\mathbf{X}    \right] 
\,
\operatorname {E}_Y \left[\mathbf{Y}  \right] 
\end{align*}
As a consequence, $\operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]$  can be written explicitly as
\begin{align*}
\operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right] 
&=
\begin{bmatrix}
\sum\limits_{i=1}^m \operatorname {E}_X \left[ X_{1i}    \right] 
\,\operatorname {E}_Y \left[  Y_i     \right] 
\\
\vdots  
\\
\sum\limits_{i=1}^m \operatorname {E}_X \left[ X_{vi}    \right] 
\,\operatorname {E}_Y \left[  Y_i     \right] 
\end{bmatrix}.
\end{align*}
The covariance matrix is then written as the expected value of the product of  $v\times 1$ vector with a $1\times v$ vector  as
\begin{align*}
&
\operatorname {E} \left[(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]   )(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]  )^{\top }\right] = 
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
\\
&\quad 
\operatorname {E} \left[
\begin{bmatrix}
\sum\limits_{i=1}^m \left(X_{1i}\,Y_{i} - \operatorname {E}_X \left[ X_{1i}    \right] 
\,\operatorname {E}_Y \left[  Y_i     \right] 
\right)
\\
\vdots  
\\
\sum\limits_{i=1}^m \left(X_{vi}\,Y_{i} - \operatorname {E}_X \left[ X_{vi}    \right] 
\,\operatorname {E}_Y \left[  Y_i     \right] 
\right)
\end{bmatrix}
\begin{bmatrix}
\sum\limits_{i=1}^m \left(X_{1i}\,Y_{i} - \operatorname {E}_X \left[ X_{1i}    \right] 
\,\operatorname {E}_Y \left[  Y_i     \right] 
\right)
\\
\vdots  
\\
\sum\limits_{i=1}^m \left(X_{vi}\,Y_{i} - \operatorname {E}_X \left[ X_{vi}    \right] 
\,\operatorname {E}_Y \left[  Y_i     \right] 
\right)
\end{bmatrix}
^{\top }\right]
.
\end{align*}
The covariance matrix has a dimension of $v\times v$
PART II
By $\operatorname {E} \left[(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]   )(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]  )^{\top }\right]_{k,l}$ I denote the element   on at the $k^\textrm{th}$ row and $l^\textrm{th}$ column of the covariance matrix. Since the expectation the covariance matrix is equal to the matrix of expecatations of the covarniance matrix' elements, and since the expectation of a sum is equal to the sum of expectations, I write $\operatorname {E} \left[(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]   )(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]  )^{\top }\right]_{k,l}$ as:
\begin{align*}
&
\operatorname {E} \left[(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]   )(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]  )^{\top }\right]_{k,l} 
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
\\
&\quad =
\sum\limits_{i=1}^m\sum\limits_{j=1}^m \operatorname {E} \left[\left(X_{ki}\,Y_{i} - \operatorname {E}_X \left[ X_{ki}    \right] 
\,\operatorname {E}_Y \left[  Y_i     \right] 
\right) \left(X_{lj}\,Y_{j} - \operatorname {E}_X \left[ X_{lj}    \right] 
\,\operatorname {E}_Y \left[  Y_j     \right] 
\right)\right]
\\
&\quad =
\sum\limits_{i=1}^m\sum\limits_{j=1}^m \operatorname {E}_X \left[ X_{ki}\,    X_{lj}\,    \right]
\,
\operatorname {E}_Y \left[ Y_{i}   \,Y_{j}   \right]
\\
&\quad -
\sum\limits_{i=1}^m\sum\limits_{j=1}^m  \operatorname {E}_X \left[ X_{ki}     \right]
\,
\operatorname {E}_X \left[ X_{lj}    \right] 
\,
\operatorname {E}_Y \left[ Y_{i}    \right]
\,
\operatorname {E}_Y \left[  Y_j     \right] 
\\
&\quad -
\sum\limits_{i=1}^m\sum\limits_{j=1}^m  \operatorname {E}_X \left[ X_{ki}    \right] 
\,
\operatorname {E}_X \left[ X_{lj}    
\right]
\,
\operatorname {E}_Y \left[  Y_i     \right]
\,
\operatorname {E}_Y \left[  Y_{j}   
\right]
\\
&\quad +
\sum\limits_{i=1}^m\sum\limits_{j=1}^m    \operatorname {E}_X \left[ X_{ki}    \right] 
\,
\operatorname {E}_X \left[ X_{lj}    \right] 
\,
\operatorname {E}_Y \left[  Y_i     \right] 
\,
\operatorname {E}_Y \left[  Y_j     \right] 
%%%%%%%%%%%55
%%%%%%%%%%%%%
%%%%%%%%%%%%
\\
&\quad =
\sum\limits_{i=1}^m\sum\limits_{j=1}^m \operatorname {E}_X \left[ X_{ki}\,    X_{lj}\,    \right]
\,
\operatorname {E}_Y \left[ Y_{i}   \,Y_{j}   \right]
\\
&\quad -
\sum\limits_{i=1}^m\sum\limits_{j=1}^m  \operatorname {E}_X \left[ X_{ki}     \right]
\,
\operatorname {E}_X \left[ X_{lj}    \right] 
\,
\operatorname {E}_Y \left[ Y_{i}    \right]
\,
\operatorname {E}_Y \left[  Y_j     \right] 
.
\end{align*}
I now attempt to separate the variables
\begin{align*}
&
\operatorname {E} \left[(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]   )(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]  )^{\top }\right]_{k,l} 
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
\\
&\quad =
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\left(
\operatorname {E}_X \left[ X_{ki}\,    X_{lj} \right]
-
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\right)
\,
\operatorname {E}_Y \left[ Y_{i}   \,Y_{j}   \right]
\\
&\quad +
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\left(
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\right)
\,
\operatorname {E}_Y \left[ Y_{i}   \,Y_{j}   \right]
\\
&\quad -
\sum\limits_{i=1}^m\sum\limits_{j=1}^m  \operatorname {E}_X \left[ X_{ki}     \right]
\,
\operatorname {E}_X \left[ X_{lj}    \right] 
\,
\operatorname {E}_Y \left[ Y_{i}    \right]
\,
\operatorname {E}_Y \left[  Y_j     \right] 
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
\\
&\quad =
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\left(
\operatorname {E}_X \left[ X_{ki}\,    X_{lj} \right]
-
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\right)
\,
\left(
\operatorname {E}_Y \left[ Y_{i}   \,Y_{j}   \right]
- 
\operatorname {E}_Y \left[ Y_{i}      \right]
\,
\operatorname {E}_Y \left[  Y_{j}   \right]
\right)
\\
&\quad +
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\left(
\operatorname {E}_X \left[ X_{ki}\,    X_{lj} \right]
-
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\right)
\,
\left( 
\operatorname {E}_Y \left[ Y_{i}      \right]
\,
\operatorname {E}_Y \left[  Y_{j}   \right]
\right)
\\
&\quad +
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\left(
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\right)
\,
\operatorname {E}_Y \left[ Y_{i}   \,Y_{j}   \right]
\\
&\quad -
\sum\limits_{i=1}^m\sum\limits_{j=1}^m  \operatorname {E}_X \left[ X_{ki}     \right]
\,
\operatorname {E}_X \left[ X_{lj}    \right] 
\,
\operatorname {E}_Y \left[ Y_{i}    \right]
\,
\operatorname {E}_Y \left[  Y_j     \right] 
.
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
\\
&\quad =
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\left(
\operatorname {E}_X \left[ X_{ki}\,    X_{lj} \right]
-
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\right)
\,
\left(
\operatorname {E}_Y \left[ Y_{i}   \,Y_{j}   \right]
- 
\operatorname {E}_Y \left[ Y_{i}      \right]
\,
\operatorname {E}_Y \left[  Y_{j}   \right]
\right)
\\
&\quad +
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\,
\operatorname {E}_Y \left[ Y_{i}      \right]
\,
\operatorname {E}_Y \left[  Y_{j}   \right]
\,
\left(
\operatorname {E}_X \left[ X_{ki}\,    X_{lj} \right]
-
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\right)
\\
&\quad +
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\,
\left(
\operatorname {E}_Y \left[ Y_{i}   \,Y_{j}   
\right]
-
\operatorname {E}_Y \left[ Y_{i}      
\right]
\,
\operatorname {E}_Y \left[ Y_{j}   
\right]
\right)
\end{align*}
From the definition of covariance [3], I rewrite the above as follows.
\begin{align*}
&
\operatorname {E} \left[(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]   )(\mathbf{X} \,\mathbf{Y}- \operatorname {E} \left[\mathbf{X} \,\mathbf{Y}  \right]  )^{\top }\right]_{k,l} 
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%
\\
&\quad =
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\operatorname {cov}_X( X_{ki},    X_{lj})
\,
\operatorname {cov}_Y( Y_{i},    Y_{j})
\\
&\quad +
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\,
\operatorname {E}_Y \left[ Y_{i}      \right]
\,
\operatorname {E}_Y \left[  Y_{j}   \right]
\,
\operatorname {cov}_X( X_{ki},    X_{lj})
\\
&\quad +
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
\,
\operatorname {cov}_Y( X_{i},    Y_{j})
%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%
\\
&\quad = 
\sum\limits_{i=1}^m\sum\limits_{j=1}^m 
\Bigl(
\operatorname {cov}_X( X_{ki},    X_{lj})
 + \operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right]
 \Bigr)\Bigl(
 \operatorname {cov}_Y( Y_{i},    Y_{j}
 )
 + \operatorname {E}_Y \left[ Y_{i}     \right]
\operatorname {E}_Y \left[ Y_{j}    \right]
 \Bigr) -\operatorname {E}_X \left[ X_{ki}     \right]
\operatorname {E}_X \left[ X_{lj}    \right] \operatorname {E}_Y \left[ Y_{i}     \right]
\operatorname {E}_Y \left[ Y_{j}    \right]
\end{align*}
BIBLIOGRAPHY
[1] https://en.wikipedia.org/wiki/Variance#For_vector-valued_random_variables 
[2] https://en.wikipedia.org/wiki/Product_distribution#Expectation_of_product_of_random_variables
[3] https://en.wikipedia.org/wiki/Covariance#Definition
 A: I would suggest amending your work by simplifying the notation and the algebra, because a clear derivation is easier to check and more convincing than a long one and puts much less burden on your readers, as well as revealing the key ideas in the result.

Let $(A,B,C,D)$ be any random variables for which $(A,C)$ is independent of $(B,D).$ From that fact and the definition of covariance as $$\operatorname{Cov}(X,Y)=E[XY]-E[X]E[Y],$$ you may compute
$$\eqalign{
\operatorname{Cov}(AB,CD) &= E[ABCD]-E[AB]E[CD] \\&= E[AC]E[BD] - E[A]E[B]E[C]E[D] \\
&= \left(\operatorname{Cov}(A,C)+E[A]E[C]\right) \left(\operatorname{Cov}(B,D)+E[B]E[D]\right) - E[A]E[B]E[C]E[D].
}$$
Consequently, letting $A=X_{ki}, B=y_i, C=X_{lj},$ and $D=y_j,$ the definition of matrix multiplication and the bilinearity of covariance yield
$$\eqalign{
\operatorname{Cov}((Xy)_k, (Xy)_l) &= \operatorname{Cov}\left(\sum_i X_{ki}y_i\ \sum_j X_{lj}y_j\right) \\
&= \sum_{i,j}\operatorname{Cov}\left(X_{ki}y_i X_{lj}y_j\right) \\
&= \sum_{i,j}\left(\operatorname{Cov}(X_{ki},X_{lj})+E[X_{ki}]E[X_{lj}]\right) \left(\operatorname{Cov}(y_i,y_j)+E[y_i]E[y_j]\right) - E[X_{ki}]E[y_i]E[X_{lj}]E[y_j],
}$$
agreeing with your result.
Because this is a straightforward calculation, few would consider it a "theorem:" formula might be a better term.
