This is an addendum to @Macro's very nice answer which lays out
exactly what needs to known in order to determine the variance of
the product of two correlated random variables. Since
\begin{align}
\operatorname{var}(XY) &= E\left[(XY)^2\right] - \left(E[XY]\right)^2
\tag{1}\\
&= E[(XY)^2] - \left(\operatorname{cov}(X,Y)+E[X]E[Y]\right)^2\\
&= E[X^2Y^2] - \left(\operatorname{cov}(X,Y)+E[X]E[Y]\right)^2\tag{2}\\
&= \left(\operatorname{cov}(X^2,Y^2)+E[X^2]E[Y^2]\right)
- \left(\operatorname{cov}(X,Y)+E[X]E[Y]\right)^2\tag{3}\\
\end{align}
where $\operatorname{cov}(X,Y)$, $E[X]$, $E[Y]$, $E[X^2]$, and
$E[Y^2]$ can be assumed to
be known quantities, we need to be able to determine the value of
$E\left[X^2Y^2\right]$ in $(2)$ or $\operatorname{cov}(X^2,Y^2)$ in $(3)$.
This is not easy to do in general, but, as pointed out already, if
$X$ and $Y$ are independent random variables, then
$\operatorname{cov}(X,Y) = \operatorname{cov}(X^2,Y^2) = 0$.
In fact, dependence, not correlation (or lack thereof) is the
key issue. That we know that $\operatorname{cov}(X,Y)$ equals $0$
instead of some nonzero value does not, by itself, help in the
least in our efforts are determining the value of
$E\left[X^2Y^2\right]$ or $\operatorname{cov}(X^2,Y^2)$ even though it
does simplify the right sides of $(2)$ and $(3)$ a little.
When $X$ and $Y$ are dependent
random variables, then in at least one (fairly common
or fairly important) special
case, it is possible to find
the value of $E\left[X^2Y^2\right]$ relatively easily.
Suppose that $X$ and $Y$ are jointly normal random variables
with correlation coefficient $\rho$. Then, conditioned
on $X = x$, the conditional density of $Y$ is a normal
density with mean
$E[Y] + \rho\left.\left.\sqrt{\frac{\operatorname{var}(Y)}{\operatorname{var}(X)}}
\right(x-E[X]\right)$ and variance $\operatorname{var}(Y)(1-\rho^2)$. Thus,
\begin{align}E[X^2Y^2 \mid X] &= X^2E[Y^2 \mid X]\\
&= X^2\left[\operatorname{var}(Y)(1-\rho^2)
+ \left(E[Y] + \rho\left.\left.\sqrt{\frac{\operatorname{var}(Y)}{\operatorname{var}(X)}}
\right(X-E[X]\right)\right)^2\right]
\end{align}
which is a quartic function of $X$, say $g(X)$, and the Law of Iterated
Expectation tells us that
$$E[X^2Y^2] = E\left[E[X^2Y^2\mid X]\right] = E[g(X)]\tag{4}$$
where the right side of $(4)$ can be computed from knowledge of the
3rd and 4th moments of $X$ -- standard results that can be found
in many texts and reference books
(meaning that I am too lazy to look them up
and include them in this answer). Fortunately, Moderator whuber has gone a step further and provided the exact result for the variance of $XY$ in a recent comment on this answer:
$$\operatorname{var}(XY) =(\mu_Y\sigma_X)^2 + (\mu_X\sigma_Y)^2 + (1+\rho^2)(\sigma_X\sigma_Y)^2 +2\rho\mu_X\mu_Y\sigma_X\sigma_Y.\tag{4}$$
Notice that when $X$ and $Y$ are independent normal random variables, Eq. $(4)$ reduces to
$$\operatorname{var}(XY) =(\mu_Y\sigma_X)^2 + (\mu_X\sigma_Y)^2 + \sigma_X^2\sigma_Y^2 \tag{5}$$
which is just what the OP and Macro wrote but expressed in slightly different notation.
Further addendum: In a now-deleted answer, @Hydrologist gives the variance of $XY$ as
$$\mathrm{Var}\left[xy\right] = \left(\mathrm{E}\left[x\right]\right)^2\mathrm{Var}\left[y\right] + \left(\mathrm{E}\left[y\right]\right)^2\mathrm{Var}\left[x\right] + 2\mathrm{E}\left[x\right]\mathrm{Cov}\left[x,y^2\right] + 2\mathrm{E}\left[y\right]\mathrm{Cov}\left[x^2,y\right]\\ + 2\mathrm{E}\left[x\right]\mathrm{E}\left[y\right]\mathrm{Cov}\left[x,y\right] +\mathrm{Cov}\left[x^2,y^2\right] - \left(\mathrm{Cov}\left[x,y\right]\right)^2 \tag{5}$$
and claims that this formula is from two papers published a half-century ago in JASA. This formula is an incorrect transcription of the results in the paper(s) cited by Hydrologist. Specifically, $\mathrm{Cov}\left[x^2,y^2\right]$ is a mistranscription of
$E[(x-E[x])^2(y-E[y])^2]$ in the journal article, and similarly for $\mathrm{Cov}\left[x^2,y\right]$ and $\mathrm{Cov}\left[x,y^2\right]$.
{\rm var}
instead of\operatorname{var}
you set a bad example. These\rm
lacks context-dependent spacing, so you'll see things like $2{\rm var}(X)$ instead of $2\operatorname{var}(X)$, and if you omit the parentheses and write $\operatorname{var}X$ then\rm
will cause you to see instead ${\rm var}X. \qquad$ $\endgroup$