从零实现时间序列奇异谱分析用Python完整复现SSA算法时间序列分析是数据科学中的核心技能之一而奇异谱分析Singular Spectrum Analysis, SSA作为一种无需预设模型的方法近年来在金融预测、气候研究、信号处理等领域展现出独特优势。本文将带您从零开始用Python完整实现SSA算法并通过可视化手段直观理解其工作原理。1. 环境准备与数据生成在开始SSA算法实现前我们需要搭建Python环境并生成一组包含趋势、周期和噪声的合成数据。这不仅能验证算法效果也为后续分析提供明确参照。import numpy as np import matplotlib.pyplot as plt plt.style.use(seaborn) # 使用更美观的绘图样式 # 参数设置 days 180 # 总天数 period 30 # 周期天数 np.random.seed(42) # 固定随机种子保证可复现性 # 生成趋势项线性下降 trend np.linspace(2, -2, numdays) # 生成周期项正弦波 time_points np.linspace(0, 2*np.pi*days/period, numdays) seasonality np.sin(time_points) * 1.5 # 生成随机噪声 noise np.random.normal(0, 0.5, days) # 组合信号 signal trend seasonality noisy_signal signal noise # 可视化各成分 fig, (ax1, ax2) plt.subplots(2, 1, figsize(12, 8)) ax1.plot(trend, labelTrend, colororange) ax1.plot(seasonality, labelSeasonality, colorgreen) ax1.plot(noise, labelNoise, alpha0.5, colorred) ax1.set_title(Signal Components) ax1.legend() ax2.plot(noisy_signal, labelNoisy Signal, colorblue) ax2.plot(signal, labelTrue Signal, linestyle--, colorblack) ax2.set_title(Composite Signal) ax2.legend() plt.tight_layout() plt.show()这段代码生成了三个关键组成部分趋势项表现为线性下降的黄色线条周期项呈现规则波动的绿色正弦曲线噪声项随机分布的红色点状波动可视化输出将清晰展示原始信号与加噪信号的对比为后续去噪效果评估提供基准。2. SSA算法核心四步实现SSA算法的核心流程可分为四个关键步骤我们将逐一实现并解释每个步骤的数学意义和实现细节。2.1 构建轨迹矩阵嵌入阶段轨迹矩阵构建是SSA的第一步其本质是将一维时间序列升维到矩阵空间。关键参数是窗口长度L的选择它直接影响分析结果。def build_trajectory_matrix(series, L): 构建轨迹矩阵 N len(series) K N - L 1 X np.zeros((L, K)) for i in range(L): X[i, :] series[i:iK] return X # 窗口长度选择通常取序列长度的1/3到1/2 L 60 traj_matrix build_trajectory_matrix(noisy_signal, L) # 可视化轨迹矩阵 plt.figure(figsize(10, 6)) plt.imshow(traj_matrix, cmapviridis, aspectauto) plt.colorbar(labelValue) plt.title(fTrajectory Matrix (L{L})) plt.xlabel(Window Position) plt.ylabel(Window Length) plt.show()窗口长度L的选择至关重要L值过小难以捕捉长周期模式L值过大可能导致过度拟合和计算效率下降经验法则通常取N/3到N/2之间其中N为序列长度2.2 SVD分解与特征分析对轨迹矩阵进行奇异值分解SVD这是SSA算法的核心数学操作能够揭示时间序列的内在结构。def perform_svd(X): 执行SVD分解并计算贡献率 U, sigma, VT np.linalg.svd(X, full_matricesFalse) # 计算各分量的能量贡献 energy sigma**2 / np.sum(sigma**2) return U, sigma, VT, energy U, sigma, VT, energy perform_svd(traj_matrix) # 绘制奇异值谱 plt.figure(figsize(12, 5)) plt.subplot(121) plt.plot(sigma, o-) plt.title(Singular Values (Log Scale)) plt.yscale(log) plt.xlabel(Component Index) plt.ylabel(Singular Value) plt.subplot(122) plt.plot(np.cumsum(energy), o-) plt.title(Cumulative Energy) plt.xlabel(Component Index) plt.ylabel(Explained Variance Ratio) plt.axhline(0.9, colorred, linestyle--, alpha0.5) plt.tight_layout() plt.show()SVD分解后我们得到三个矩阵U矩阵代表时间域的特征模式Σ矩阵对角矩阵包含奇异值反映各成分重要性V矩阵代表序列域的特征模式奇异值谱图可帮助我们确定重要成分的数量通常前几个成分对应趋势和主要周期。2.3 成分分组与重构根据奇异值谱的分析结果我们需要对成分进行合理分组这是SSA最具技巧性的环节。def reconstruct_component(U, sigma, VT, component_idx, L): 重构指定成分 s np.zeros(len(sigma)) s[component_idx] sigma[component_idx] S np.diag(s) X_comp U S VT return diagonal_averaging(X_comp, L) def diagonal_averaging(X, L): 对角平均将矩阵转换回时间序列 N X.shape[1] L - 1 reconstructed np.zeros(N) for k in range(N): if k L - 1: reconstructed[k] np.mean([X[i, k-i] for i in range(k1)]) elif k X.shape[1]: reconstructed[k] np.mean([X[i, k-i] for i in range(L)]) else: reconstructed[k] np.mean([X[i, k-i] for i in range(k-X.shape[1]1, L)]) return reconstructed # 重构前4个主要成分 components [] for i in range(4): components.append(reconstruct_component(U, sigma, VT, i, L)) # 可视化各成分 plt.figure(figsize(12, 8)) for i, comp in enumerate(components): plt.subplot(4, 1, i1) plt.plot(comp) plt.title(fComponent {i1} (Energy: {energy[i]:.2%})) plt.tight_layout() plt.show()分组策略通常包括趋势成分通常对应第一个奇异值周期成分对应后续几个显著奇异值噪声成分剩余的小奇异值2.4 结果分析与可视化最后一步是将重构的成分与原始信号进行对比评估SSA的分解效果。# 组合趋势和周期成分 reconstructed_trend components[0] reconstructed_seasonality components[1] components[2] # 可视化对比 plt.figure(figsize(12, 8)) plt.subplot(311) plt.plot(trend, labelTrue Trend, colorblack, linestyle--) plt.plot(reconstructed_trend, labelReconstructed Trend) plt.title(Trend Component Comparison) plt.legend() plt.subplot(312) plt.plot(seasonality, labelTrue Seasonality, colorblack, linestyle--) plt.plot(reconstructed_seasonality, labelReconstructed Seasonality) plt.title(Seasonal Component Comparison) plt.legend() plt.subplot(313) plt.plot(signal, labelTrue Signal, colorblack, linestyle--) plt.plot(reconstructed_trend reconstructed_seasonality, labelReconstructed Signal) plt.title(Full Signal Comparison) plt.legend() plt.tight_layout() plt.show()通过对比图可以直观评估趋势成分的捕捉是否准确周期信号的振幅和相位是否匹配噪声抑制的效果如何3. 参数优化与实用技巧SSA算法的效果很大程度上依赖于参数选择和操作技巧本节将分享实际应用中的经验。3.1 窗口长度L的优化选择窗口长度L是SSA最重要的参数没有统一选择标准但有一些实用准则选择方法说明适用场景1/3规则取N/3N为序列长度通用场景周期倍数取主要周期的整数倍已知周期特征试错法尝试多个L值比较结果复杂序列# 测试不同L值的效果 L_values [30, 60, 90, 120] fig, axes plt.subplots(len(L_values), 1, figsize(12, 10)) for ax, L in zip(axes, L_values): traj_matrix build_trajectory_matrix(noisy_signal, L) U, sigma, VT, _ perform_svd(traj_matrix) comp1 reconstruct_component(U, sigma, VT, 0, L) ax.plot(trend, labelTrue Trend, linestyle--) ax.plot(comp1, labelfL{L}) ax.legend() ax.set_title(fWindow Length L {L}) plt.tight_layout() plt.show()实验表明L60序列长度的1/3在本案例中表现最佳既能捕捉趋势又不会引入过多噪声。3.2 成分选择与噪声处理如何选择有效成分、舍弃噪声成分是SSA应用的关键决策点。常用方法包括能量阈值法保留累计能量达到90%的成分奇异值拐点法选择奇异值下降曲线的拐点频率分析法通过FFT检查各成分的频率特征# 成分能量分析示例 cum_energy np.cumsum(energy) threshold 0.9 # 90%能量阈值 n_components np.argmax(cum_energy threshold) 1 print(f保留前{n_components}个成分可解释{threshold:.0%}的方差)在实际项目中建议结合多种方法综合判断并通过业务知识验证成分的合理性。4. 完整SSA算法封装与扩展为了便于实际应用我们将前述步骤封装为完整的Python类并添加实用功能。class SingularSpectrumAnalysis: def __init__(self, series, window_lengthNone): self.series series self.N len(series) self.L window_length if window_length else self.N // 3 self.K self.N - self.L 1 def decompose(self, n_componentsNone): 执行完整SSA分解 # 构建轨迹矩阵 self.X build_trajectory_matrix(self.series, self.L) # SVD分解 self.U, self.sigma, self.VT, self.energy perform_svd(self.X) # 确定成分数量 if n_components is None: cum_energy np.cumsum(self.energy) n_components np.argmax(cum_energy 0.9) 1 # 重构成分 self.components [] for i in range(n_components): comp reconstruct_component(self.U, self.sigma, self.VT, i, self.L) self.components.append(comp) return self.components def reconstruct(self, component_indices): 重构指定成分组合 return np.sum([self.components[i] for i in component_indices], axis0) def plot_components(self, n_components6): 可视化主要成分 plt.figure(figsize(12, 8)) for i in range(min(n_components, len(self.components))): plt.subplot(n_components, 1, i1) plt.plot(self.components[i]) plt.title(fComponent {i1} (Energy: {self.energy[i]:.2%})) plt.tight_layout() plt.show() # 使用示例 ssa SingularSpectrumAnalysis(noisy_signal, L60) components ssa.decompose() ssa.plot_components() # 重构信号 reconstructed ssa.reconstruct([0, 1, 2]) # 组合趋势和前两个周期成分 # 性能评估 from sklearn.metrics import mean_squared_error mse mean_squared_error(signal, reconstructed) print(fReconstruction MSE: {mse:.4f})这个封装类提供了以下扩展功能自动窗口长度选择智能成分数量确定灵活的成分组合重构内置可视化方法重构质量评估指标在实际项目中可以进一步扩展以下功能自动化成分分组算法实时更新机制处理流数据结合预测模型进行时间序列预测并行计算加速大规模数据处理
保姆级教程：用Python的NumPy和Matplotlib一步步拆解时间序列（含SSA算法完整代码）
从零实现时间序列奇异谱分析用Python完整复现SSA算法时间序列分析是数据科学中的核心技能之一而奇异谱分析Singular Spectrum Analysis, SSA作为一种无需预设模型的方法近年来在金融预测、气候研究、信号处理等领域展现出独特优势。本文将带您从零开始用Python完整实现SSA算法并通过可视化手段直观理解其工作原理。1. 环境准备与数据生成在开始SSA算法实现前我们需要搭建Python环境并生成一组包含趋势、周期和噪声的合成数据。这不仅能验证算法效果也为后续分析提供明确参照。import numpy as np import matplotlib.pyplot as plt plt.style.use(seaborn) # 使用更美观的绘图样式 # 参数设置 days 180 # 总天数 period 30 # 周期天数 np.random.seed(42) # 固定随机种子保证可复现性 # 生成趋势项线性下降 trend np.linspace(2, -2, numdays) # 生成周期项正弦波 time_points np.linspace(0, 2*np.pi*days/period, numdays) seasonality np.sin(time_points) * 1.5 # 生成随机噪声 noise np.random.normal(0, 0.5, days) # 组合信号 signal trend seasonality noisy_signal signal noise # 可视化各成分 fig, (ax1, ax2) plt.subplots(2, 1, figsize(12, 8)) ax1.plot(trend, labelTrend, colororange) ax1.plot(seasonality, labelSeasonality, colorgreen) ax1.plot(noise, labelNoise, alpha0.5, colorred) ax1.set_title(Signal Components) ax1.legend() ax2.plot(noisy_signal, labelNoisy Signal, colorblue) ax2.plot(signal, labelTrue Signal, linestyle--, colorblack) ax2.set_title(Composite Signal) ax2.legend() plt.tight_layout() plt.show()这段代码生成了三个关键组成部分趋势项表现为线性下降的黄色线条周期项呈现规则波动的绿色正弦曲线噪声项随机分布的红色点状波动可视化输出将清晰展示原始信号与加噪信号的对比为后续去噪效果评估提供基准。2. SSA算法核心四步实现SSA算法的核心流程可分为四个关键步骤我们将逐一实现并解释每个步骤的数学意义和实现细节。2.1 构建轨迹矩阵嵌入阶段轨迹矩阵构建是SSA的第一步其本质是将一维时间序列升维到矩阵空间。关键参数是窗口长度L的选择它直接影响分析结果。def build_trajectory_matrix(series, L): 构建轨迹矩阵 N len(series) K N - L 1 X np.zeros((L, K)) for i in range(L): X[i, :] series[i:iK] return X # 窗口长度选择通常取序列长度的1/3到1/2 L 60 traj_matrix build_trajectory_matrix(noisy_signal, L) # 可视化轨迹矩阵 plt.figure(figsize(10, 6)) plt.imshow(traj_matrix, cmapviridis, aspectauto) plt.colorbar(labelValue) plt.title(fTrajectory Matrix (L{L})) plt.xlabel(Window Position) plt.ylabel(Window Length) plt.show()窗口长度L的选择至关重要L值过小难以捕捉长周期模式L值过大可能导致过度拟合和计算效率下降经验法则通常取N/3到N/2之间其中N为序列长度2.2 SVD分解与特征分析对轨迹矩阵进行奇异值分解SVD这是SSA算法的核心数学操作能够揭示时间序列的内在结构。def perform_svd(X): 执行SVD分解并计算贡献率 U, sigma, VT np.linalg.svd(X, full_matricesFalse) # 计算各分量的能量贡献 energy sigma**2 / np.sum(sigma**2) return U, sigma, VT, energy U, sigma, VT, energy perform_svd(traj_matrix) # 绘制奇异值谱 plt.figure(figsize(12, 5)) plt.subplot(121) plt.plot(sigma, o-) plt.title(Singular Values (Log Scale)) plt.yscale(log) plt.xlabel(Component Index) plt.ylabel(Singular Value) plt.subplot(122) plt.plot(np.cumsum(energy), o-) plt.title(Cumulative Energy) plt.xlabel(Component Index) plt.ylabel(Explained Variance Ratio) plt.axhline(0.9, colorred, linestyle--, alpha0.5) plt.tight_layout() plt.show()SVD分解后我们得到三个矩阵U矩阵代表时间域的特征模式Σ矩阵对角矩阵包含奇异值反映各成分重要性V矩阵代表序列域的特征模式奇异值谱图可帮助我们确定重要成分的数量通常前几个成分对应趋势和主要周期。2.3 成分分组与重构根据奇异值谱的分析结果我们需要对成分进行合理分组这是SSA最具技巧性的环节。def reconstruct_component(U, sigma, VT, component_idx, L): 重构指定成分 s np.zeros(len(sigma)) s[component_idx] sigma[component_idx] S np.diag(s) X_comp U S VT return diagonal_averaging(X_comp, L) def diagonal_averaging(X, L): 对角平均将矩阵转换回时间序列 N X.shape[1] L - 1 reconstructed np.zeros(N) for k in range(N): if k L - 1: reconstructed[k] np.mean([X[i, k-i] for i in range(k1)]) elif k X.shape[1]: reconstructed[k] np.mean([X[i, k-i] for i in range(L)]) else: reconstructed[k] np.mean([X[i, k-i] for i in range(k-X.shape[1]1, L)]) return reconstructed # 重构前4个主要成分 components [] for i in range(4): components.append(reconstruct_component(U, sigma, VT, i, L)) # 可视化各成分 plt.figure(figsize(12, 8)) for i, comp in enumerate(components): plt.subplot(4, 1, i1) plt.plot(comp) plt.title(fComponent {i1} (Energy: {energy[i]:.2%})) plt.tight_layout() plt.show()分组策略通常包括趋势成分通常对应第一个奇异值周期成分对应后续几个显著奇异值噪声成分剩余的小奇异值2.4 结果分析与可视化最后一步是将重构的成分与原始信号进行对比评估SSA的分解效果。# 组合趋势和周期成分 reconstructed_trend components[0] reconstructed_seasonality components[1] components[2] # 可视化对比 plt.figure(figsize(12, 8)) plt.subplot(311) plt.plot(trend, labelTrue Trend, colorblack, linestyle--) plt.plot(reconstructed_trend, labelReconstructed Trend) plt.title(Trend Component Comparison) plt.legend() plt.subplot(312) plt.plot(seasonality, labelTrue Seasonality, colorblack, linestyle--) plt.plot(reconstructed_seasonality, labelReconstructed Seasonality) plt.title(Seasonal Component Comparison) plt.legend() plt.subplot(313) plt.plot(signal, labelTrue Signal, colorblack, linestyle--) plt.plot(reconstructed_trend reconstructed_seasonality, labelReconstructed Signal) plt.title(Full Signal Comparison) plt.legend() plt.tight_layout() plt.show()通过对比图可以直观评估趋势成分的捕捉是否准确周期信号的振幅和相位是否匹配噪声抑制的效果如何3. 参数优化与实用技巧SSA算法的效果很大程度上依赖于参数选择和操作技巧本节将分享实际应用中的经验。3.1 窗口长度L的优化选择窗口长度L是SSA最重要的参数没有统一选择标准但有一些实用准则选择方法说明适用场景1/3规则取N/3N为序列长度通用场景周期倍数取主要周期的整数倍已知周期特征试错法尝试多个L值比较结果复杂序列# 测试不同L值的效果 L_values [30, 60, 90, 120] fig, axes plt.subplots(len(L_values), 1, figsize(12, 10)) for ax, L in zip(axes, L_values): traj_matrix build_trajectory_matrix(noisy_signal, L) U, sigma, VT, _ perform_svd(traj_matrix) comp1 reconstruct_component(U, sigma, VT, 0, L) ax.plot(trend, labelTrue Trend, linestyle--) ax.plot(comp1, labelfL{L}) ax.legend() ax.set_title(fWindow Length L {L}) plt.tight_layout() plt.show()实验表明L60序列长度的1/3在本案例中表现最佳既能捕捉趋势又不会引入过多噪声。3.2 成分选择与噪声处理如何选择有效成分、舍弃噪声成分是SSA应用的关键决策点。常用方法包括能量阈值法保留累计能量达到90%的成分奇异值拐点法选择奇异值下降曲线的拐点频率分析法通过FFT检查各成分的频率特征# 成分能量分析示例 cum_energy np.cumsum(energy) threshold 0.9 # 90%能量阈值 n_components np.argmax(cum_energy threshold) 1 print(f保留前{n_components}个成分可解释{threshold:.0%}的方差)在实际项目中建议结合多种方法综合判断并通过业务知识验证成分的合理性。4. 完整SSA算法封装与扩展为了便于实际应用我们将前述步骤封装为完整的Python类并添加实用功能。class SingularSpectrumAnalysis: def __init__(self, series, window_lengthNone): self.series series self.N len(series) self.L window_length if window_length else self.N // 3 self.K self.N - self.L 1 def decompose(self, n_componentsNone): 执行完整SSA分解 # 构建轨迹矩阵 self.X build_trajectory_matrix(self.series, self.L) # SVD分解 self.U, self.sigma, self.VT, self.energy perform_svd(self.X) # 确定成分数量 if n_components is None: cum_energy np.cumsum(self.energy) n_components np.argmax(cum_energy 0.9) 1 # 重构成分 self.components [] for i in range(n_components): comp reconstruct_component(self.U, self.sigma, self.VT, i, self.L) self.components.append(comp) return self.components def reconstruct(self, component_indices): 重构指定成分组合 return np.sum([self.components[i] for i in component_indices], axis0) def plot_components(self, n_components6): 可视化主要成分 plt.figure(figsize(12, 8)) for i in range(min(n_components, len(self.components))): plt.subplot(n_components, 1, i1) plt.plot(self.components[i]) plt.title(fComponent {i1} (Energy: {self.energy[i]:.2%})) plt.tight_layout() plt.show() # 使用示例 ssa SingularSpectrumAnalysis(noisy_signal, L60) components ssa.decompose() ssa.plot_components() # 重构信号 reconstructed ssa.reconstruct([0, 1, 2]) # 组合趋势和前两个周期成分 # 性能评估 from sklearn.metrics import mean_squared_error mse mean_squared_error(signal, reconstructed) print(fReconstruction MSE: {mse:.4f})这个封装类提供了以下扩展功能自动窗口长度选择智能成分数量确定灵活的成分组合重构内置可视化方法重构质量评估指标在实际项目中可以进一步扩展以下功能自动化成分分组算法实时更新机制处理流数据结合预测模型进行时间序列预测并行计算加速大规模数据处理
