Power spectral density of the output of a linear time invariant system with a weakly-stationary process Please note that I have no background whatsoever in Fourier analysis and very little in time series.
This is exercise 2.10 in Theodoridis' Machine Learning:

Show that the autocorrelation of the output of a linear system, with impulse response $w_{n}$, $n \in \mathbb{Z}$, is related to the autocorrelation of the intput WSS process, via
  $$r_d(k) = r_u(k) * w_k * w^{*}_{-k}$$

I don't even know where to get started with this exercise. Here is the information I've gathered:
Some notation/terms to define:


*

*WSS = Wide-Sense Stationary (or weakly stationary)

*$r_d(k) = \mathbb{E}[d_nd_{n-k}]$

*$r_u(k) = \mathbb{E}[u_nu_{n-k}]$

*$w^{*}$ means the complex conjugate of $w$ (I think)

*$*$ refers to the convolution sum. For example, if $\dots, w_0, w_1, \dots$ are parameters and $(u_n)$ a time series, $$w_n * u_n = \sum_{i=-\infty}^{+\infty}w^{*}_i u_{n-i}$$ I'm not sure how this works when you have two convolutions (are convolutions associative?).


It is worth noting that this exercise is used as an intermediate step to the proof of theorem 

Theorem 2.1. The power spectral density of the output, $d_n$, of a linear time invariant system, where it is excited by a WSS stochastic process, $u_n$, is given by
  $$S_d(\omega) = |W(\omega)|^2S_u(\omega)$$

The proof is solved by using the result of the exercise above and then, apparently, using the "well-known" properties of the Fourier transform implying that
$$r_u(k) * w_k \mapsto S_u(\omega)W(\omega)$$
$$w^{*}_{-k} \mapsto W^{*}(\omega)$$
I opened one of my Time Series texts and found this theorem in the text by Shumway and Stoffer (4th ed.) on p. 175 (with very different notation), but I don't understand the proof, given my lack of time series and Fourier analysis background.
Could someone please assist with helping me start exercise 2.10, addressing my concerns in bold above?
 A: 
Show that $$r_d = r_u \star w \star \tilde w \tag{1}$$ where $r_d(k)$ is the autocorrelation function of the WSS output process, $r_u(k)$ is the autocorrelation function of the WSS input process, $w_k$ is the impulse response of the linear time-invariant system, $\tilde w_k = w_{-k}$ and
  $\star$ denotes the convolution operator.

This is just a matter of massaging the sums involved and identifying the various components in the result.  Define the discrete-time autocorrelation function of the deterministic sequence $w$ as 
$$r_w(k) = \sum_{n=-\infty}^\infty w_nw_{n-k} \tag{2}$$ which, since 
$w_{n-k} = \tilde{w}_{k-n}$, can be expressed in terms of $\star$ as
$$r_w(k) = \sum_{n=-\infty}^\infty w_n\tilde w_{k-n} = w\star \tilde w$$ so that the result to be proven can be expressed as
$$r_d = r_u \star (w \star \tilde w) = r_u \star r_w. \tag{3}$$
Note that since $r_w$ is an even function -- we have that $r_w(k) = r_w(-k)$ for all $k$ -- the convolution $(3)$ can be expressed as
$$r_d(k) = \sum_{n=-\infty}^\infty r_u(n)r_w(k-n) = \sum_{n=-\infty}^\infty r_u(n)r_w(n-k), \tag{4}$$
that is, $r_d$ is also the cross-correlation of the functions $r_u$ and $r_w$.
Now, note that
$\displaystyle d_n = \sum_{j=-\infty}^\infty u_j w_{n-j}$
and so
\begin{align} r_d(k) &= E[d_nd_{n-k}]\\
&= E\left[\sum_{j=-\infty}^\infty u_j w_{n-j}\sum_{m=-\infty}^\infty u_m w_{n-k-m}\right]\\
&= \sum_{j=-\infty}^\infty \sum_{m=-\infty}^\infty E[u_ju_m]  w_{n-j}w_{n-k-m}\\
&= \sum_{j=-\infty}^\infty \sum_{m=-\infty}^\infty r_u(j-m) w_{n-j}w_{n-k-m}
\end{align}
The standard way of thinking of this double sum over $j,m \in \mathbb Z$ is that we first compute the inner sum (over $m$) for a fixed value of $j$, and then compute the outer sum (on $j$), that is, we work on column sums first and then on the rows of this infinite array. We get the same result if we fix $s = j-m$ and sum over all $j,m$ such that $j-m=s$, and then vary $s$; that is, we work on a diagonal of this array in computing the "inner" sum, and then sum over all possible diagonals. Thus we have
\begin{align}
r_d(k) &= \sum_{j=-\infty}^\infty \sum_{m=-\infty}^\infty r_u(j-m) w_{n-j}w_{n-k-m}\\
&= \sum_{s=-\infty}^\infty r_u(s)\left[ \sum_{j,m\colon j-m=s}^\infty  w_{n-j}w_{n-k-m}\right]\\
&= \sum_{s=-\infty}^\infty r_u(s) \left[\sum_{j=-\infty}^\infty  w_{n-j}w_{(n-j)+j-k-m}\right]&\scriptstyle{\text{Now set}~t=n-j, j-m=s}\\
&= \sum_{s=-\infty}^\infty r_u(s) \left[\sum_{t=-\infty}^\infty  w_{t}w_{t-k+s}\right]\\
&= \sum_{s=-\infty}^\infty r_u(s) r_w(k-s)\\
\text{that is,}\qquad r_d &= r_u \star r_w,
\end{align}
thus proving $(3)$.
A: You are very close to the solution, basically you are missing one step.
Your linear system is defined by
$$ d(n) = u(n)*w(n) $$
Where $u(n)$ is a WSS. This means that the convolutions is how you model the effect of going through the linear system modeled by $w(n)$. 
Now consider that the Fourier transform of a convolution in time is equivalent to a multiplication in frequency space.
$$ D(f) = U(f)W(f) $$
Being the $D(f)$ Fourier transform of the output
$$ D(f) = F \{u(n)\} $$
Consider as well that the Fourier transform of an auto-correlation function is the power spectral density.
Now, apply this principle to your equations above
$$ r_d(k) = r_u(k)*w(k)*w^*(-k) $$
$$ F\{r_d(k)\} = S_d(f) = S_u(f) W(f) W^*(f) $$
$$ S_d(f) = S_u(f) |W(f)|^2 $$ 
Reaching Theorem 2.1 and proving the equality right.
As for your question  are convolutions associative , I guess you meant linear: yes they are.
$$ (u_1 + u_2) * w  = u_1*w + u_2*w $$
Think of it as integrals, they are linear.
