Let $(X_n)$ and $(Y_n)$ be independent wide-sense stationary (WSS) processes with averages $\mu_X=EX_n$ and $\mu_Y=EY_n$, covariance functions $C_X(m)=E(X_{n+m}-EX_{n+m})(X_{n}-EX_{n})$ and $C_Y(m)=E(Y_{n+m}-EY_{n+m})(Y_{n}-EY_{n})$, correlation functions $R_X(m)=C_X(m)/C_X(0)$ and $R_Y(m)=C_Y(m)/C_Y(0)$, and power spectral densities (PSDs) $f$ and $g$, so that
$R_X(m)=\int_{-\infty}^\infty e^{2\pi i mx} f(x)dx$ and $R_Y(m)=\int_{-\infty}^\infty e^{2\pi i mx} g(x)dx$. Let $Z_n=X_n Y_n$.

Then the process $(Z_n)$ is also WSS, with the average
$$\mu_Z=EX_n Y_n=EX_n\,EY_n=\mu_X \mu_Y,$$
covariance function
\begin{align*}
C_Z(m)&=E(X_{n+m}Y_{n+m}-\mu_X\mu_Y)(X_{n}Y_{n}-\mu_X\mu_Y) \\
&=EX_{n+m}X_{n}\,EY_{n+m}Y_n-\mu_X^2\mu_Y^2\\
&=(C_X(m)+\mu_X^2)(C_Y(m)+\mu_Y^2)-\mu_X^2\mu_Y^2\\
&=C_X(m)C_Y(m)+\mu_Y^2C_X(m)+\mu_X^2C_Y(m),
\end{align*}
correlation function
\begin{equation*}
R_Z(m)=\frac{C_Z(m)}{C_Z(0)}
=a\,R_X(m)R_Y(m)
+b\,R_X(m)
+c\,R_Y(m)
\end{equation*}
with $a:=\frac{C_X(0)C_Y(0)}{C_Z(0)}$, $b:=\frac{C_X(0)\mu_Y^2}{C_Z(0)}$, $c:=\frac{C_Y(0)\mu_X^2}{C_Z(0)}$,
and PSD
\begin{equation*}
h=a f*g+bf+cg,
\end{equation*}
where $f*g$ is the convolution of $f$ and $g$.