difference in the spectral densities of autoregressive sequence by FFT and analytical solution

Question

I want to obtain the power spectral density (PSD) of an autoregressive sequence, AR(1). The analytical solution according to this reference (page 12) is

For $X_t = \phi_1X_{t-1}+W_t, W_t \sim N(0,\sigma^2_w)$，

$$f(w) = \frac{\sigma_w^2}{1-2\phi_1 \cos(2 \pi w)+\phi_1^2}$$

Let $\sigma_w = 1$, $\phi_1 = 0.4$, then the corresponding PSD curve in [0, 2pi] using fft and analytical solution is

As can be seen from the figure below, the PSD obtained from fft is far from being a U-shaped curve obtained by the analytical solution.

Can anyone enlighten me on why the results are different?

---

The code I used to generate the plot is:

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq  

def analy(x):
    PSD = 1/(1.16-0.8*np.cos(2*np.pi*x))
    return PSD

def my_fft(N):   
    T = 1  # sample spacing
    t = np.linspace(0.0, N*T, N, endpoint=True)
    y = np.zeros([N,1])

    rho = 0.4
    for i in range(len(t)-1):
        y[i+1] = rho*y[i] + (1-0)*np.random.normal(0,1) 

    yf = fft(y)   
    xf = fftfreq(N, T)[:N//2]  # sample frequency points
    PSD = 2.0/N * np.abs(yf[0:N//2])

    return xf*np.pi*2, PSD

def compare_PSD():
    N = 100     # Number of sample points
    UB = np.pi  # sample upper bound
    x_analy = np.linspace(0, UB, num=N)
    PSD_analy = analy(x_analy)
    x_fft, PSD_fft = my_fft(N)
    plt.plot(x_analy, PSD_analy, label='Analytical')
    plt.plot(x_fft, PSD_fft, label='scipy.fft')
    plt.xlabel('Frequency')
    plt.ylabel('PSD')
    plt.legend()
    plt.show()

compare_PSD() 

Answer 1

I got the correct results. There are two mistakes in my previous comparison.

1. Conceptual error: PSD is obtained by taking the Fourier transformation on the covariance, not the original time series data.

$$ \begin{aligned} f(v) &=\sum_{h=-\infty}^{\infty} \frac{\sigma_{w}^{2}}{1-\phi_{1}^{2}} \cdot \rho(h) e^{-2 \pi i v h} \\ &=\frac{\sigma_{w}^{2}}{1-\phi_{1}^{2}}\left[2 \sum_{h \in(0, \infty)} \phi_{1}^{h} e^{-2 \pi i v h}- \phi_1^{0} e^{0}\right] \end{aligned} $$

2. Coding: "scipy.fft" takes a vector as its input. However, "y" in the previous code is a matrix.

Below is the correct result.

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft, fftfreq  

def cal_PSD_anal(x, phi):
    PSD = 1/(1 + phi**2-2*phi*np.cos(2*np.pi*x))
    return PSD

def cal_PSD_fft(N, phi):     
    # covaraince 
    gamma = phi**np.arange(N)
    yf = fft(gamma).real
    PSD = 1 / (1 - phi**2) * (2*np.abs(yf[0:N//2]) - 1)    
    return PSD

def compare_PSD():
    N = 100     # Number of sample points
    UB = 0.5 #np.pi  # sample upper bound
    T = 1  # sample spacing
    phi = 0.4
        
    x = fftfreq(N, T)[:N//2]  # sample frequency points
    PSD_fft = cal_PSD_fft(N, phi)    

    PSD_anal = cal_PSD_anal(x, phi)
    
    plt.figure(figsize=(4, 3), dpi=300) 
    plt.plot(x, PSD_anal, 'bo', label='Analytical', markersize=4)
    plt.plot(x, PSD_fft, 'r', label='FFT') 

    plt.title('PSD of AR(1, {})'.format(phi))
    plt.xlabel('Frequency')
    plt.ylabel('PSD')
    plt.legend()
    plt.grid(axis='y')
    plt.tight_layout()
    plt.savefig('./PSD_AR_1_{}.png'.format(phi), dpi=300)
    plt.show()
    
compare_PSD()