季节性与趋势建模：傅里叶谐波与样条拟合

y(t) = a0 + Σ [aₙ cos(2π n t / P) + bₙ sin(2π n t / P)]

其中 `n` 为谐波阶数，阶数越高，能捕获的季节形态越精细。

### 如何在回归中使用

将序列的时间索引 `t`（如1,2,3…）或日期特征转换为若干对正弦、余弦变量，作为回归模型的输入。对于长度为 `P` 的季节周期，第 `k` 阶谐波对应两个特征：

- `cos(2π k t / P)`
- `sin(2π k t / P)`

通常取 `K` 个谐波（推荐 `K = 2` 或 `3` 对于月度季度数据，`K = 3` 到 `6` 对于周数据）就能很好地近似常见季节形状。

### Python 示例：生成傅里叶特征

```python
import numpy as np
import pandas as pd

def fourier_features(t, period, n_harmonics):
    """
    t: 时间索引 (1D array)
    period: 季节周期长度
    n_harmonics: 保留的谐波对数
    返回 DataFrame，列名包含 cos/sin 与阶数
    """
    features = {}
    t = np.array(t)
    for k in range(1, n_harmonics + 1):
        features[f'cos_{k}'] = np.cos(2 * np.pi * k * t / period)
        features[f'sin_{k}'] = np.sin(2 * np.pi * k * t / period)
    return pd.DataFrame(features)

# 使用示例：模拟 3 年的月度数据（周期 12）
dates = pd.date_range('2020-01-01', periods=36, freq='M')
t = np.arange(1, 37)
harmonics = fourier_features(t, period=12, n_harmonics=3)
print(harmonics.head())

选择谐波数量

• 过少：季节形状被过度平滑，丢失波峰波谷细节。
• 过多：可能拟合噪声，导致过拟合。
• 经验建议：先用 2～3 个谐波，观察残差是否仍存在周期性模式；若有，可逐步增加。交叉验证也能辅助选择。

样条函数与趋势建模

样条（Spline）是一种分段多项式，在节点（knots）处满足连续性条件。常用的是 三次样条，它保证函数值、一阶导和二阶导在节点处连续，因此曲线非常光滑。

回归样条

在回归模型中，可将时间变量通过样条基函数展开，例如自然三次样条的基函数 B(t)，然后最小二乘拟合：

trend(t) = Σ βⱼ Bⱼ(t)

节点位置控制趋势的灵活性：节点越多，曲线越能弯曲，但也越容易过拟合。

Python 实现：使用 patsy 或 sklearn 的样条基

推荐使用 patsy 的内置函数，或 scipy.interpolate.BSpline，或者简单的自然样条变换来自定义。这里展示用 patsy 的 cr() 函数（需要安装 patsy）：

import patsy
import pandas as pd

# 假设 df 包含时间列 t
# 自然三次样条，内部节点在 25% 和 75% 分位数
transformed = patsy.dmatrix("cr(t, df=3)", {"t": t}, return_type='dataframe')

df 参数控制自由度，即最终基函数的个数。df=3 通常给出一个很平滑的趋势（仅两个内部节点）。你可以增大 df 来增加灵活性。

若不想依赖 patsy，也可以手动生成 B 样条基，例如使用 scipy.interpolate.BSpline 并构建设计矩阵。

节点选择策略

• 等距节点：适用于时间序列无明显结构变化。
• 分位数节点：在数据密集处放更多节点，处理变化剧烈区域。
• 手动定节点：根据业务里程碑（如政策变化、产品发布）设定。
• 正则化光滑样条（如 Ridge 惩罚）自动控制平滑度，避免手动调 df。

集成模型：趋势 + 季节 + 残差

将趋势成分与季节成分放入一个加法模型（或乘法模型经对数变换后变为加法）：

y(t) = trend(t) + seasonality(t) + error

其中趋势由样条基表示，季节由傅里叶谐波表示。这实际上就是一个线性回归模型，可以很方便地用 OLS（普通最小二乘）拟合。

完整建模步骤

1. 准备时间特征：生成整数时间索引 t。
2. 生成傅里叶特征：确定周期 P（如年=12，周=7），选择谐波数 K。
3. 生成样条基：选择样条类型（自然三次样条）和自由度 df。
4. 合并设计矩阵：包含截距项（通常自动包含）、所有样条基列、所有正余弦列。
5. 拟合线性模型：使用 statsmodels 或 sklearn.linear_model.LinearRegression。
6. 诊断残差：检查残差自相关性（如 ACF 图），必要时调整谐波数或样条自由度。
7. 预测：将未来的 t 带入模型，提前计算好对应的样条基和傅里叶特征，得到预测值。

完整代码示例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression

# 模拟数据：趋势 + 年季节 + 噪声
np.random.seed(42)
dates = pd.date_range('2018-01-01', periods=60, freq='M')
t = np.arange(1, 61)
trend = 0.05 * t + 0.0005 * t**2  # 轻微二次趋势
seasonal = 2 * np.sin(2*np.pi*t/12) + 1.5 * np.cos(2*np.pi*t/12)
y = trend + seasonal + np.random.normal(0, 0.8, size=len(t))

df = pd.DataFrame({'y': y, 't': t})

# 1. 傅里叶特征 (周期12, 3个谐波)
def add_fourier(df, period, n_harmonics, col_prefix=''):
    t = df['t'].values
    for k in range(1, n_harmonics+1):
        df[f'{col_prefix}cos_{k}'] = np.cos(2*np.pi*k*t/period)
        df[f'{col_prefix}sin_{k}'] = np.sin(2*np.pi*k*t/period)
    return df

df = add_fourier(df, period=12, n_harmonics=3)

# 2. 样条基 (自然三次样条, df=4)
# 使用 patsy 生成 (需要安装: pip install patsy)
import patsy
spline_basis = patsy.dmatrix("cr(t, df=4)", {"t": df['t']}, return_type='dataframe')
# 为避免共线性，dmatrix 自动包含截距，我们稍后去掉自己的截距
# 将样条列（去掉第一列截距）加入 df
spline_cols = [f'spline_{i}' for i in range(spline_basis.shape[1]-1)]
for i, col in enumerate(spline_cols):
    df[col] = spline_basis.iloc[:, i+1]

# 3. 设计矩阵 X (不包含额外截距，因为 LinearRegression 会自动加)
feature_cols = spline_cols + [c for c in df.columns if 'cos_' in c or 'sin_' in c]
X = df[feature_cols]
y = df['y']

# 4. 拟合模型
model = LinearRegression()
model.fit(X, y)
print(f"R² score: {model.score(X, y):.4f}")

# 5. 查看拟合与分解
df['trend_hat'] = model.predict(X.assign(**{c: 0 for c in df.columns if c in feature_cols and 'cos' in c or 'sin' in c}))
# 上面代码不直接，更清晰做法：单独预测趋势（季节特征置0）
X_trend = X.copy()
for c in X_trend.columns:
    if 'cos' in c or 'sin' in c:
        X_trend[c] = 0.0
df['trend_hat'] = model.predict(X_trend)

df['season_hat'] = model.predict(X) - df['trend_hat']

# 绘图
plt.figure(figsize=(12,6))
plt.plot(df['t'], y, label='Original')
plt.plot(df['t'], model.predict(X), label='Fitted', linewidth=2)
plt.plot(df['t'], df['trend_hat'], '--', label='Trend')
plt.legend()
plt.title('Trend + Seasonal Fit')
plt.show()