Let's consider independent random vectors $\hat{\boldsymbol\theta}_i$, $i = 1, \dots, m$, which are all unbiased for $\boldsymbol\theta$ and that $$\mathbb{E}\left[\left(\hat{\boldsymbol\theta}_i - \boldsymbol\theta\right)^{T}\left(\hat{\boldsymbol\theta}_i - \boldsymbol\theta\right)\right] = \sigma^2\text{.}$$ Let $\mathbf{1}_{n \times p}$ be the $n \times p$ matrix of all ones.
Consider the problem of finding $$\mathbb{E}\left[\left(\hat{\boldsymbol\theta} - \boldsymbol\theta\right)^{T}\left(\hat{\boldsymbol\theta} - \boldsymbol\theta\right)\right]$$ where $$\hat{\boldsymbol\theta} = \dfrac{1}{m}\sum_{i=1}^{m}\hat{\boldsymbol\theta}_i\text{.}$$
My attempt is to notice the fact that $$\hat{\boldsymbol\theta} = \dfrac{1}{m}\underbrace{\begin{bmatrix} \hat{\boldsymbol\theta}_1 & \hat{\boldsymbol\theta}_2 & \cdots & \hat{\boldsymbol\theta}_m \end{bmatrix}}_{\mathbf{S}}\mathbf{1}_{m \times 1}$$ and thus $$\text{Var}(\hat{\boldsymbol\theta}) = \dfrac{1}{m^2}\text{Var}(\mathbf{S}\mathbf{1}_{m \times 1})\text{.}$$ How does one find the variance of a random matrix times a constant vector? You may assume that I am familiar with finding variances of linear transformations of a random vector: i.e., if $\mathbf{x}$ is a random vector, $\mathbf{b}$ a vector of constants, and $\mathbf{A}$ a matrix of constants, assuming all are comformable, $$\mathbb{E}[\mathbf{A}\mathbf{x}+\mathbf{b}] = \mathbf{A}\mathbb{E}[\mathbf{x}]+\mathbf{b}$$ $$\mathrm{Var}\left(\mathbf{A}\mathbf{x}+\mathbf{b}\right)=\mathbf{A}\mathrm{Var}(\mathbf{x})\mathbf{A}^{\prime}$$