A general approach is to directly parameterize the covariance matrix, then fit a Gausssian distribution using this parameterization (e.g. using maximum likelihood or a Bayesian method).
For example, suppose the covariance matrix $C$ (of size $d \times d$) is constrained to have a constant variance $\sigma^2$ along the diagonal and constant covarariance $c$ on the off-diagonal elements. This can be written as $C = c \mathbf{1} + (\sigma^2-c) I$, where $\mathbf{1}$ is a matrix of ones and $I$ is the identity matrix. We require that $\sigma^2 \ge c$, guaranteeing that $C$ is positive semidefinite.
I found maximum likelihood estimators for the variance and covariance in this case, which turn out to be simple averages. The variance $\hat{\sigma}^2$ is obtained by averaging the sample variance for each variable. And, the covariance $\hat{c}$ is obtained by averaging the sample covariance for each pair of non-identical variables. That is, suppose the ordinary (unconstrained) sample covariance matrix is $\tilde{C} = \frac{1}{n} \sum_{i=1}^n (x_i-\mu) (x_i-\mu)^T$. Then:
$$\hat{\sigma}^2 = \frac{1}{d} \tilde{C}_{ii} \quad \quad
\hat{c} = \frac{1}{d (d-1)} \sum_{i=1}^d \sum_{j \ne i} \tilde{C}_{ij}$$
Derivation
Likelihood function
Assuming the data $X = \{x_1, \dots, x_n\} \subset \mathbb{R}^d$ are i.i.d. Gaussian, the likelihood given mean $\mu$ and covariance matrix $C$ is:
$$p(X \mid \mu, C) = \prod_{i=1}^n
(2 \pi)^{-\frac{d}{2}} \det(C)^{-\frac{1}{2}}
\exp \left( -\frac{1}{2} (x_i-\mu)^T C^{-1} (x_i-\mu) \right)$$
The special form of $C$ above allows some simplifications. In particular, its inverse and determinant are:
$$C^{-1} = a \mathbf{1} + b I \quad \quad
\det(C) = \frac{1}{b^d + d a b^{d-1}}$$
$$\text{where:} \quad
a = \frac{c}{(c - \sigma^2) (d c - c + \sigma^2)} \quad \quad
b = \frac{1}{\sigma^2 - c}$$
After a little algebra, we can write the negative log likelihood $\mathcal{L}$ as a function of the parameters:
$$\mathcal{L}(\mu, \sigma^2, c) =
\frac{n d}{2} \log (2 \pi)
- \frac{n}{2} \log(b^d + d a b^{d-1})
+ \frac{a}{2} S_1 + \frac{b}{2} S_2$$
$$\text{where:} \quad
S_1 = \sum_{i=1}^n (x_i-\mu)^T \mathbf{1} (x_i-\mu) \quad \quad
S_2 = \sum_{i=1}^n (x_i-\mu)^T (x_i-\mu)$$
Maximum likelihood estimation
Maximum likelihood parameter estimates are obtained by minimizing $\mathcal{L}$. Closed form solutions can be found by differentiating $\mathcal{L}$ with respect to $\sigma^2$ and $c$, setting the derivatives to zero, and solving. This yields:
$$\hat{\sigma}^2 = \frac{1}{n d} \sum_{i=1}^n (x_i-\mu)^T (x_i-\mu)
\quad \quad
\hat{c} = \frac{1}{n d (d-1)}
\sum_{i=1}^n \left[
(x_i-\mu)^T (\mathbf{1} - I) (x_i-\mu)
\right]$$
The maximum likelihood estimate for $\mu$ is just the sample mean of the data, and can be plugged into $\hat{\sigma}^2$ and $\hat{c}$.
Writing things out in scalar form makes it apparent that the variance $\hat{\sigma}^2$ simply averages the sample variance of all variables. And, the covariance $\hat{c}$ simply averages the sample covariance between each pair of non-identical variables.
