text:  covered$$ f(\omega) = \gamma(0) + 2 \sum_{k \geq 1} \gamma(k) \cos(\omega k), \qquad \omega \in \mathbb
  $$ I(\omega) := \frac{1}{n} \left| \sum_{t=0}^{n-1} X_t e^{i t \omega} \right|^2, \qquad \omega \in \mathbb R \tag{1} $$
  (Recall that $ |z| $ denotes the modulus of complex number $ z $)
  Alternatively, $ I(\omega) $ can be expressed as$$ I(\omega) = \frac{1}{n} \left\{ \left[\sum_{t=0}^{n-1} X_t \cos(\omega t) \right]^2 + \left[\sum_{t=0}^{n-1} X_t \sin(\omega t) \right]^2 \right\} $$$$ \omega_j := \frac{2 \pi j}{n}, \quad j = 0, \ldots, n - 1 $$
  In what sense is $ I(\omega_j) $ an estimate of $ f(\omega_j) $?
  The answer is straightforward, although it does involve some algebra.
  With a bit of effort, one can show that for any integer $ j > 0 $,$$ \sum_{t=0}^{n-1} e^{i t \omega_j } = \sum_{t=0}^{n-1} \exp \left\{ i 2 \pi j \frac{t}{n} \right\} = 0 $$
  Letting $ \bar X $ denote the sample mean $ n^{-1} \sum_{t=0}^{n-1} X_t $, we then have$$ n I(\omega_j) = \left| \sum_{t=0}^{n-1} (X_t - \bar X) e^{i t \omega_j } \right|^2 = \sum_{t=0}^{n-1} (X_t - \bar X) e^{i t \omega_j } \sum_{r=0}^{n-1} (X_r - \bar X) e^{-i r \omega_j } $$
  By carefully working through the sums, one can transform this to$$ n I(\omega_j) = \sum_{t=0}^{n-1} (X_t - \bar X)^2 + 2 \sum_{k=1}^{n-1} \sum_{t=k}^{n-1} (X_t - \bar X)(X_{t-k} - \bar X) \cos(\omega_j k) $$
  Now let$$ \hat \gamma(k) := \frac{1}{n} \sum_{t=k}^{n-1} (X_t - \bar X)(X_{t-k} - \bar X), \qquad k = 0,1,\ldots, n-1 $$$$ I(\omega_j) = \hat \gamma(0) + 2 \sum_{k=1}^{n-1} \hat \gamma(k) \cos(\omega_j k) $$
  statsmodels.tsa.stattools.periodogram in the
  statsmodels$$ A_j := \sum_{t=0}^{n-1} a_t \exp \left\{ i 2 \pi \frac{tj}{n} \right\}, \qquad j = 0, \ldots, n-1 $$
  With
  numpy.fft.fft imported as
  fft and $ a_0, \ldots, a_{n-1} $ stored in NumPy array
  a, the function call
  fft(a) returns the values $ A_0, \ldots, A_{n-1} $ as a NumPy array.
  It follows that when the data $ X_0, \ldots, X_{n-1} $ are stored in array
  X, the values $ I(\omega_j) $ at the Fourier frequencies, which are given by
  can be computed by
  np.abs(fft(X))**2 / len(X).
  Note: The NumPy function
  abs acts elementwise, and correctly handles complex numbers (by computing their modulus, which is exactly what we need).
  A function called
  periodogram that puts all this together can be found here.
  Let’s generate some data for this function using the
  ARMA class from QuantEcon.py (see the lecture on linear processes for more details).
  Here’s a code snippet that, once the preceding code has been run, generates data from the process
  $$ X_t = 0.5 X_{t-1} + \epsilon_t - 0.8 \epsilon_{t-2} \tag{2} $$
  where $ \{ \epsilon_t \} $ is white noise with unit variance, and compares the periodogram to the actual spectral density
  import matplotlib.pyplot as plt %matplotlib inline from quantecon import ARMA, periodogram n = 40 # Data size ϕ, θ = 0.5, (0, -0.8) # AR and MA parameters lp = ARMA(ϕ, θ) X = lp.simulation(ts_length=n) fig, ax = plt.subplots() x, y = periodogram(X) ax.plot(x, y, 'b-', lw=2, alpha=0.5, label='periodogram') x_sd, y_sd = lp.spectral_density(two_pi=False, res=120) ax.plot(x_sd, y_sd, 'r-', lw=2, alpha=0.8, label='spectral density') ax.legend() plt.show()
  /home/qebuild/anaconda3/lib/python3.7/site-packages/numpy/core/numeric.py:538: ComplexWarning: Casting complex values to real discards the imaginary part return array(a, dtype, copy=False, order=order) parametric assumptions.
  Typically, nonparametric estimation of densities requires some degree of smoothing.
  The standard way that smoothing is applied to periodograms is by taking local averages.
  In other words, the value $ I(\omega_j) $ is replaced with a weighted average of the adjacent values$$ I(\omega_{j-p}), I(\omega_{j-p+1}), \ldots, I(\omega_j), \ldots, I(\omega_{j+p}) $$
  This weighted average can be written as
  $$ I_S(\omega_j) := \sum_{\ell = -p}^{p} w(\ell) I(\omega_{j+\ell}) \tag{3} $$
  where the weights $ w(-p), \ldots, w(p) $ are a sequence of $ 2p + 1 $ nonnegative values summing to one.
  In general,).
  import numpy as np def hanning_window(M): w = [0.5 - 0.5 * np.cos(2 * np.pi * n/(M-1)) for n in range(M)] return w window = hanning_window(25) / np.abs(sum(hanning_window(25))) x = np.linspace(-12, 12, 25) fig, ax = plt.subplots(figsize=(9, 7)) ax.plot(x, window) ax.set_title("Hanning window") ax.set_ylabel("Weights") ax.set_xlabel("Position in sequence of weights") plt.show()
  Estimation with Smoothing¶
  Our next step is to provide code that will not only estimate the periodogram but also provide smoothing as required.
  Such functions have been written in estspec.py and are available once you’ve installed QuantEcon.py.
  np.convolve.
  Readers are left either to explore or simply to use this code according to their interests.
  The next three figures each show smoothed and unsmoothed periodograms, as well as the population or “true” spectral density.
  (The model is the same as before — see equation (2) — and there are 400 observations)
  From.py. are assumed to be generated by an AR(1) process.
  (More general ARMA settings can be handled using similar techniques to those described below)
  Suppose in particular that $ \{X_t\} $ is covariance stationary and AR(1), with
  $$ X_{t+1} = \mu + \phi X_t + \epsilon_{t+1} \tag{4} $$
  $$ f(\omega) = \left| \frac{1}{1 - \phi e^{i\omega}} \right|^2 g(\omega) \tag{5} $$
  This suggests that the recoloring step, which constructs an estimate $ I $ of $ f $ from $ I_0 $, should set$$ I(\omega) = \left| \frac{1}{1 - \hat \phi e^{i\omega}} \right|^2 I_0(\omega) $$
  where $ \hat \phi $ is the OLS estimate of $ \phi $.
  The code for
  ar_periodogram() — the third function in
  estspec.py —.
  Exercise 1¶
  Replicate this figure (modulo randomness).
  The model is as in equation (2) and there are 400 observations.
  For the smoothed periodogram, the window.
  from quantecon import ar_periodogram
  ## Data n = 400 ϕ = 0.5 θ = 0, -0.8 lp = ARMA(ϕ, θ) X = lp.simulation(ts_length=n) fig, ax = plt.subplots(3, 1, figsize=(10, 12)) for i, wl in enumerate((15, 55, 175)): # window lengths x, y = periodogram(X) ax[i].plot(x, y, 'b-', lw=2, alpha=0.5, label='periodogram') x_sd, y_sd = lp.spectral_density(two_pi=False, res=120) ax[i].plot(x_sd, y_sd, 'r-', lw=2, alpha=0.8, label='spectral density') x, y_smoothed = periodogram(X, window='hamming', window_len=wl) ax[i].plot(x, y_smoothed, 'k-', lw=2, label='smoothed periodogram') ax[i].legend() ax[i].set_title(f'window length = {wl}') plt.show()
  lp = ARMA(-0.9) wl = 65 fig, ax = plt.subplots(3, 1, figsize=(10,12)) for i in range(3): X = lp.simulation(ts_length=150) ax[i].set_xlim(0, np.pi) x_sd, y_sd = lp.spectral_density(two_pi=False, res=180) ax[i].semilogy(x_sd, y_sd, 'r-', lw=2, alpha=0.75, label='spectral density') x, y_smoothed = periodogram(X, window='hamming', window_len=wl) ax[i].semilogy(x, y_smoothed, 'k-', lw=2, alpha=0.75, label='standard smoothed periodogram') x, y_ar = ar_periodogram(X, window='hamming', window_len=wl) ax[i].semilogy(x, y_ar, 'b-', lw=2, alpha=0.75, label='AR smoothed periodogram') ax[i].legend(loc='upper left') plt.show()