时间序列系数聚类流程

1. 肘部法选 K → 2. K-Means 聚类 → 3. 3D 可视化
2. 简单说：把长得像、变化规律像的时间序列分到一类。

3. 常用两种路子：

   (1). 基于原始数据直接聚类（不推荐）

   直接拿整条序列算距离，容易被趋势、漂移、幅度干扰，结果不准。

   (2). 基于特征提取聚类

   步骤：（1）对每条时间序列拟合模型（AR、ARIMA 等）（2）提取系数 / 特征（提取 AR或ARMA 系数）（3）对这些系数向量做 K-Means 聚类（4）系数相近 → 序列动态特征相近 → 算同一类

4. 案例：12个美国工业生产指数（非耐用品制造业），月度指标

   library('forecast')
   # install.packages("rgl")
   library('rgl')   #画3D图
   data=read.csv("data.csv",header=TRUE)
   data=data[3:467,13:24]
   datanew=data[2:465,]
   colnames(datanew)<-data[1,]
   datanew=as.matrix(datanew)

   max_lag=15   # 设定AR模型的最大滞后阶数p=15
   coef=matrix(0,16,12)  # 新建16行12列的0矩阵，存12个序列的AR系数
   for(j in 1:12)  # 循环：逐列处理（共12列数据）
   {
     # 核心：拟合 ARIMA(p,1,0) 模型 = 一阶差分后的 AR(p)
     fit=auto.arima(as.numeric(datanew[,j]),  # 第j列数据转数值
                    max.p=max_lag, # p最大=15
                    max.q=0,d=1,  # 强制q=0，不做MA
                    ic='bic',  # 用BIC信息准则选阶数
                    allowdrift = FALSE)  # **不允许漂移项（drift）
     print(fit$coef) # 打印拟合出的系数
     # 如果系数非空，就存入矩阵
     if(length(fit$coef)!=0)
     {
       coef[(1:length(fit$coef)),j]=fit$coef
     }
   }
   coef=coef[1:3,]   # 只保留**前3行**（前3个AR系数）
   coef=t(coef)  # 矩阵转置 → 最终 12行3列

   
   #####k-means#####
   # 1. 初始化一个长度为10的向量，用来存每个K对应的SSE（组内误差平方和）
   SSE=numeric(10)
   # 2. 循环：尝试 K=1 到 K=10 个聚类
   for(k in 1:10)
   {
     # 对 coef 矩阵做 K-Means 聚类
     # nstart=20：随机初始20次，选最优结果（避免局部最优）
     fit_kmeans=kmeans(x=coef,k,nstart = 20)
     # 把【组内误差平方和】存进 SSE 向量
     # tot.withinss = 所有点到自己簇中心的距离平方和（越小越好）
     SSE[k]=fit_kmeans$tot.withinss
   }
   # 3. 画肘部法则图：看 SSE 随 K 增加的下降曲线
   plot(1:10,SSE,
        xlab='聚类数目',
        ylab='误差平方和SSE',
        type='b',  # 点+线都画
        main='K均值聚类图')
   abline(v=5,lty=2)

   # 4. 正式用 K=5 做最终聚类
   #nstart=20：多次随机初始，保证聚类结果稳定
   fit_kmeans=kmeans(x=coef,5,nstart = 20)
   # 输出聚类结果：簇中心、每个样本属于哪一类、SSE等
   fit_kmeans
   # 5. 3D 画聚类结果（coef 是 3列数据 → 3维图）
   # 颜色按聚类分组 +1 是为了颜色从1开始，更美观
   plot3d(coef, col = fit_kmeans$cluster+1,
          size=2,   # 点大小
          type='s',  # 画球形点
          xlab='',ylab='',zlab='')  # 不显示轴标签