首页课程笔记统计系应用时间序列分析

# 频谱密度

# 频谱密度的定义与性质

  • 频谱密度函数 f(λ)f(\lambda)：对于一个平稳时间序列的自协方差函数 γh\gamma_h，其频谱密度定义为：

    f(λ)=12π∑n=−∞∞γneinλ, λ∈[−π,π]f(\lambda)=\frac{1}{2\pi}\sum_{n=-\infty}^\infty \gamma_n e^{in\lambda},~\lambda\in [-\pi,\pi]

    它描述了时间序列的方差（或总能量）在不同频率上的分布情况。

  • 频谱密度与频谱分布函数：频谱密度是频谱分布函数 F(λ)F(\lambda) 的导数，二者关系如下：

    F(λ)=∫−πλf(ν)dνF(\lambda)=\int_{-\pi}^\lambda f(\nu)d\nu

  • 频谱分布函数与自协方差函数：自协方差函数 γh\gamma_h 可以通过频谱分布函数求得，这称为谱表示定理：

    γh=∫−ππeihνdF(ν), h=0,±1,±2,⋯\gamma_h=\int_{-\pi}^\pi e^{ih\nu}dF(\nu),~h=0,\pm1,\pm2,\cdots

# 线性时间序列的频谱密度

如果一个弱平稳、零均值时间序列 {Xt}\{X_t\} 是由另一个弱平稳、零均值时间序列 {Yt}\{Y_t\} 经过线性滤波得到的，即：

Xt=∑j=−∞∞ψjYt−j, 且 ∑j=−∞∞∣ψj∣<∞X_t=\sum_{j=-\infty}^\infty\psi_jY_{t-j},~\text{且}~\sum_{j=-\infty}^\infty|\psi_j|<\infty

则 {Xt}\{X_t\} 的频谱密度 fX(λ)f_X(\lambda) 与 {Yt}\{Y_t\} 的频谱密度 fY(λ)f_Y(\lambda) 存在如下关系：

fX(λ)=∣ψ(e−iλ)∣2fY(λ)f_X(\lambda)=|\psi(e^{-i\lambda})|^2f_Y(\lambda)

其中，ψ(B)=∑j=−∞∞ψjBj\psi(B)=\sum_{j=-\infty}^\infty \psi_j B^j 是滤波器的转移函数。

# ARMA 过程的频谱密度

对于一个 ARMA(p,q)ARMA(p,q) 过程 {yt}\{y_t\}，其方程为：

ϕ(B)yt=θ(B)εt, εt∼WN(0,σ2)\phi(B)y_t=\theta(B)\varepsilon_t,~\varepsilon_t\sim WN(0,\sigma^2)

其中 ϕ(z)=1−ϕ1z−⋯−ϕpzp\phi(z)=1-\phi_1z-\cdots-\phi_pz^p 和 θ(z)=1+θ1z+⋯+θqzq\theta(z)=1+\theta_1z+\cdots+\theta_qz^q 且 ϕpθq≠0\phi_p\theta_q\ne0。若 ϕ(z)\phi(z) 在单位圆上没有零点，则该过程的频谱密度为：

fy(λ)=σ22π∣θ(e−iλ)∣2∣ϕ(e−iλ)∣2, −π≤λ≤πf_y(\lambda)=\frac{\sigma^2}{2\pi}\frac{|\theta(e^{-i\lambda})|^2}{|\phi(e^{-i\lambda})|^2},~-\pi\le\lambda\le\pi

频谱密度的近似：对于一个实序列 {yt}\{y_t\}，其连续频谱 ff 可以被因果的 AR(p)AR(p) 过程或可逆的 MA(q)MA(q) 过程的频谱密度所近似。这表明 AR(p)AR(p) 或 MA(q)MA(q) 模型可以用来近似描述 {yt}\{y_t\}。

# 周期图估计

# 周期图的定义

周期图 In(λ)I_n(\lambda) 是一种谱密度估计方法，定义为：

In(λ)=12πn∣∑k=1nxke−ikλ∣2I_n(\lambda)=\frac{1}{2\pi n}\left|\sum_{k=1}^nx_ke^{-ik\lambda}\right|^2

其与基于样本自协方差函数（SACF）的谱密度估计 f^(λ)\hat{f}(\lambda) 有紧密联系。

基于 SACF 的谱密度估计 f^(λ)\hat{f}(\lambda)：

f^(λ)=12π∑k=−(n−1)n−1γ^ke−ikλ\hat{f}(\lambda)=\frac{1}{2\pi}\sum_{k=-(n-1)}^{n-1}\hat{\gamma}_ke^{-ik\lambda}

其中，样本自协方差函数 γ^k\hat{\gamma}_k 为：

γ^±k=1n∑t=1n−k(xt−xˉ)(xt+k−xˉ)\hat{\gamma}_{\pm k}=\frac{1}{n}\sum_{t=1}^{n-k}(x_t-\bar{x})(x_{t+k}-\bar{x})

如果将序列中心化，即 xkx_k 替换为 xk−xˉx_k-\bar{x}，则 f^(λ)\hat{f}(\lambda) 与周期图 In(λ)I_n(\lambda) 在频率 λj=2πj/n\lambda_j=2\pi j/n 处相等：

In(λj)=12πn∣∑k=1n(xk−xˉ)e−ikλj∣2=Jn(λj)I_n(\lambda_j)=\frac{1}{2\pi n}\left|\sum_{k=1}^n(x_k-\bar{x})e^{-ik\lambda_j}\right|^2 = J_n(\lambda_j)

# 周期图的性质

  • 周期图与 SACF 的关系：样本自协方差函数 γ^k\hat{\gamma}_k 可以通过对周期图进行傅里叶变换得到：

    γ^k=∫−ππIn(s)eiksds,    ∣k∣≤n−1\hat{\gamma}_k=\int_{-\pi}^\pi I_n(s)e^{iks}ds,~~~~|k|\le n-1

  • 周期图的收敛性：对于零均值弱平稳序列 {xt}\{x_t\}，如果 ∑k=1∞∣γk∣<∞\sum_{k=1}^\infty|\gamma_k|<\infty，则周期图是频谱密度的渐近无偏估计：

    E[In(λ)]→f(λ), n→∞E[I_n(\lambda)]\to f(\lambda),~n\to\infty

    但周期图的估计通常不稳定，因为其方差不收敛到零：

    Var(In(λ))≈f(λ)2≠0Var(I_n(\lambda))\approx f(\lambda)^2 \ne 0

    这主要是由于求和项过多，且滞后阶数 kk 接近 nn 时，γ^k\hat{\gamma}_k 的估计偏差增大。

# 谱密度平滑估计

为了解决周期图方差不稳定的问题，通常使用时间窗口（或滞后窗）对样本自协方差函数进行平滑，以降低高频噪声和估计误差的积累。

  • 滞后窗估计：

    f^(λ)=12π∑k=−(n−1)n−1λn(k)γ^ke−ikλ\hat{f}(\lambda)=\frac{1}{2\pi}\sum_{k=-(n-1)}^{n-1}\lambda_n(k)\hat{\gamma}_ke^{-ik\lambda}

    其中，{λn(k)}\{\lambda_n(k)\} 是时间窗口函数。对于 ARMA 模型，一种常见的时间窗口是：

    λn(k)={1∣k∣≤n0∣k∣>n\lambda_n(k)=\begin{cases}1&|k|\le\sqrt{n}\\0&|k|>\sqrt{n}\end{cases}

  • 频谱窗口 Wn(λ)W_n(\lambda)：是时间窗口 λn(k)\lambda_n(k) 的傅里叶变换：

    Wn(λ)=12π∑k=−(n−1)n−1λn(k)e−ikλW_n(\lambda)=\frac{1}{2\pi}\sum_{k=-(n-1)}^{n-1}\lambda_n(k)e^{-ik\lambda}

    它满足 λn(k)=∫−ππWn(s)eiksds\lambda_n(k)=\int_{-\pi}^\pi W_n(s)e^{iks}ds。

  • 平滑估计的卷积形式：平滑后的谱密度估计可以看作是周期图与频谱窗口的卷积：

    f^(λ)=∫−ππIn(s)Wn(λ−s)ds=In(λ)∗Wn(λ)\hat{f}(\lambda)=\int_{-\pi}^\pi I_n(s) W_n(\lambda-s)ds=I_n(\lambda)*W_n(\lambda)

    为了使估计有效，频谱窗口通常需要满足以下性质：∫−ππWn(λ)dλ=1\int_{-\pi}^\pi W_n(\lambda)d\lambda=1（即 λ0=1\lambda_0=1），∫−ππWn2(λ)dλ<∞\int_{-\pi}^\pi W^2_n(\lambda)d\lambda<\infty，Wn(−λ)=Wn(λ)W_n(-\lambda)=W_n(\lambda) 等。

# 常用窗口函数

  • 矩形窗口：

    λn(k)={1∣k∣≤mn0∣k∣>mn\lambda_n(k)=\begin{cases}1&|k|\le m_n\\0&|k|>m_n\end{cases}

    其中 mn=A[n], 1≤A≤3m_n=A[\sqrt{n}],~1\le A\le 3。其频谱窗口为Dirichlet 核：

    Wn(λ)=12πsin⁡((2mn+1)λ/2)sin⁡(λ/2)W_n(\lambda)=\frac{1}{2\pi}\frac{\sin((2m_n+1)\lambda/2)}{\sin(\lambda/2)}

  • Bartlett 窗口：

    λn(k)={1−∣k∣/mn∣k∣≤mn0∣k∣>mn\lambda_n(k)=\begin{cases}1-|k|/m_n&|k|\le m_n\\0&|k|>m_n\end{cases}

    其频谱窗口为Fejer 核：

    Wn(λ)=12πmn(sin⁡(mnλ/2)sin⁡(λ/2))2W_n(\lambda)=\frac{1}{2\pi m_n}\left(\frac{\sin(m_n\lambda/2)}{\sin(\lambda/2)}\right)^2