[proofplan]
The proof reduces the periodogram to the squared modulus of the limiting complex normal Fourier ordinate. The asymptotic normality theorem for Fourier ordinates gives convergence of $d_n(\omega_n)$ to a complex Gaussian random variable with independent real and imaginary parts. Applying the continuous mapping principle to the squared-modulus map identifies the limit as one half of a chi-squared random variable with two degrees of freedom, hence as an exponential random variable with mean $1$. Since the limiting distribution is non-degenerate, convergence in probability to the deterministic value $f(\omega)$ is impossible.
[/proofplan]
[step:Use Fourier ordinate asymptotic normality at the nearby frequencies]
By the assumed asymptotic normality theorem for Fourier ordinates, applied at the sequence of Fourier frequencies $(\omega_n)_{n \in \mathbb{N}}$ with $\omega_n \to \omega \in (0,\pi)$, there exist real-valued random variables $U$ and $V$ such that
\begin{align*}
d_n(\omega_n) \xrightarrow{d} Z := U+iV,
\end{align*}
where $U$ and $V$ are independent Gaussian random variables with distributions
\begin{align*}
U \sim \mathcal{N}\left(0,\frac{f(\omega)}{2}\right),
\qquad
V \sim \mathcal{N}\left(0,\frac{f(\omega)}{2}\right).
\end{align*}
The hypotheses required for that theorem are precisely the stationarity and regularity assumptions imposed in the statement, together with the fact that $\omega$ is an interior frequency and that $f$ is continuous at $\omega$.
[guided]
The periodogram is defined as the squared modulus of the Fourier ordinate, so the distributional behaviour of $I_n(\omega_n)$ is inherited from the distributional behaviour of $d_n(\omega_n)$. The preceding Fourier-ordinate central limit theorem applies because $\omega \in (0,\pi)$ is an interior frequency, the frequencies $\omega_n = 2\pi j_n/n$ satisfy $\omega_n \to \omega$, and the spectral density $f$ is continuous at $\omega$.
Thus there are real-valued random variables $U$ and $V$ such that
\begin{align*}
d_n(\omega_n) \xrightarrow{d} Z := U+iV,
\end{align*}
where $U$ and $V$ are independent and
\begin{align*}
U \sim \mathcal{N}\left(0,\frac{f(\omega)}{2}\right),
\qquad
V \sim \mathcal{N}\left(0,\frac{f(\omega)}{2}\right).
\end{align*}
The positivity assumption $f(\omega)>0$ matters here because it ensures that the limiting Gaussian distribution has nonzero variance in both real and imaginary directions.
[/guided]
[/step]
[step:Pass from the Fourier ordinate to the periodogram by the squared-modulus map]
Define the continuous map
\begin{align*}
\Phi: \mathbb{C} &\to [0,\infty) \\
z &\mapsto \frac{|z|^2}{f(\omega)}.
\end{align*}
Since $f(\omega)>0$, $\Phi$ is well-defined and continuous. By the continuous mapping principle (citing a result not yet in the wiki: Continuous Mapping Theorem), the convergence $d_n(\omega_n) \xrightarrow{d} Z$ implies
\begin{align*}
\frac{I_n(\omega_n)}{f(\omega)}
=
\Phi(d_n(\omega_n))
\xrightarrow{d}
\Phi(Z)
=
\frac{|Z|^2}{f(\omega)}
=
\frac{U^2+V^2}{f(\omega)}.
\end{align*}
[guided]
The periodogram is exactly the squared modulus:
\begin{align*}
I_n(\omega_n) = |d_n(\omega_n)|^2.
\end{align*}
Therefore we introduce the map
\begin{align*}
\Phi: \mathbb{C} &\to [0,\infty) \\
z &\mapsto \frac{|z|^2}{f(\omega)}.
\end{align*}
This map is continuous because $z \mapsto |z|^2$ is continuous on $\mathbb{C}$, and division by $f(\omega)$ is valid because $f(\omega)>0$.
We now apply the continuous mapping principle (citing a result not yet in the wiki: Continuous Mapping Theorem). Its hypothesis is convergence in distribution of the input random variables, which we have from
\begin{align*}
d_n(\omega_n) \xrightarrow{d} Z.
\end{align*}
Its conclusion is convergence in distribution after applying the continuous map $\Phi$:
\begin{align*}
\Phi(d_n(\omega_n)) \xrightarrow{d} \Phi(Z).
\end{align*}
Substituting the definitions of $\Phi$, $I_n$, and $Z=U+iV$ gives
\begin{align*}
\frac{I_n(\omega_n)}{f(\omega)}
=
\Phi(d_n(\omega_n))
\xrightarrow{d}
\Phi(Z)
=
\frac{|Z|^2}{f(\omega)}
=
\frac{U^2+V^2}{f(\omega)}.
\end{align*}
[/guided]
[/step]
[step:Identify the squared complex Gaussian limit as an exponential random variable]
Define standardized real-valued random variables
\begin{align*}
A := \frac{U}{\sqrt{f(\omega)/2}},
\qquad
B := \frac{V}{\sqrt{f(\omega)/2}}.
\end{align*}
Since $U$ and $V$ are independent Gaussian random variables with variance $f(\omega)/2$, the random variables $A$ and $B$ are independent and each has distribution $\mathcal{N}(0,1)$. Hence
\begin{align*}
A^2+B^2 \sim \chi^2_2.
\end{align*}
Moreover,
\begin{align*}
\frac{U^2+V^2}{f(\omega)}
=
\frac{1}{2}\left(A^2+B^2\right).
\end{align*}
The distribution $\frac{1}{2}\chi^2_2$ is $\operatorname{Exp}(1)$, so
\begin{align*}
\frac{I_n(\omega_n)}{f(\omega)} \xrightarrow{d} \operatorname{Exp}(1).
\end{align*}
[guided]
To identify the limiting distribution, we standardize the two real Gaussian components. Define
\begin{align*}
A := \frac{U}{\sqrt{f(\omega)/2}},
\qquad
B := \frac{V}{\sqrt{f(\omega)/2}}.
\end{align*}
This is well-defined because $f(\omega)>0$. Since
\begin{align*}
U \sim \mathcal{N}\left(0,\frac{f(\omega)}{2}\right),
\qquad
V \sim \mathcal{N}\left(0,\frac{f(\omega)}{2}\right),
\end{align*}
the standardized variables satisfy
\begin{align*}
A \sim \mathcal{N}(0,1),
\qquad
B \sim \mathcal{N}(0,1).
\end{align*}
Because $U$ and $V$ are independent, the rescaled variables $A$ and $B$ are also independent.
By the definition of the chi-squared distribution with two degrees of freedom, the sum of squares of two independent standard normal random variables has distribution $\chi^2_2$:
\begin{align*}
A^2+B^2 \sim \chi^2_2.
\end{align*}
Now rewrite the limiting squared modulus in terms of $A$ and $B$:
\begin{align*}
U^2 = \frac{f(\omega)}{2}A^2,
\qquad
V^2 = \frac{f(\omega)}{2}B^2,
\end{align*}
and therefore
\begin{align*}
\frac{U^2+V^2}{f(\omega)}
=
\frac{1}{2}\left(A^2+B^2\right).
\end{align*}
A chi-squared random variable with two degrees of freedom has the same distribution as $2Y$, where $Y \sim \operatorname{Exp}(1)$. Equivalently,
\begin{align*}
\frac{1}{2}\chi^2_2 \sim \operatorname{Exp}(1).
\end{align*}
Combining this identification with the convergence from the previous step yields
\begin{align*}
\frac{I_n(\omega_n)}{f(\omega)} \xrightarrow{d} \operatorname{Exp}(1).
\end{align*}
[/guided]
[/step]
[step:Exclude convergence in probability to the spectral density]
Suppose, for contradiction, that
\begin{align*}
I_n(\omega_n) \xrightarrow{\mathbb{P}} f(\omega).
\end{align*}
Since $f(\omega)>0$, division by $f(\omega)$ gives
\begin{align*}
\frac{I_n(\omega_n)}{f(\omega)} \xrightarrow{\mathbb{P}} 1.
\end{align*}
Convergence in probability to a constant implies convergence in distribution to that constant, so
\begin{align*}
\frac{I_n(\omega_n)}{f(\omega)} \xrightarrow{d} 1.
\end{align*}
But the previous step gives convergence in distribution to $\operatorname{Exp}(1)$, and the distribution $\operatorname{Exp}(1)$ is not the point mass at $1$; for example,
\begin{align*}
\mathbb{P}(Y \leq 1) = 1-e^{-1} \neq 1
\end{align*}
when $Y \sim \operatorname{Exp}(1)$. This contradicts uniqueness of limits in distribution. Therefore $I_n(\omega_n)$ does not converge in probability to $f(\omega)$.
[/step]
