R语言时间序列学习+代码记录1：趋势与季节性

数据来源主要为FRED（课程要求使用），代码部分来自于老师，参考用书：

1.Diebold FX Forecasting in Economics, Business, Finance and Beyond. University of Pennsylvania, 2017.
        - available on-line http:/​​/​​www.ssc.upenn.edu/​​~fdiebold/​​Teaching221/​​Forecasting.pdf
2.Diebold FX Elements of Forecasting, 4th edition (Cengage Publishing, 2007).
        - hard copy can be purchased used from on-line book stores
        - available on-line http:/​​/​​www.ssc.upenn.edu/​​~fdiebold/​​Teaching221/​​BookPhotocopy.pdf
        - author's version http://www.ssc.upenn .edu/​​~fdiebold/​​Teaching221/​​FullBook.pdf

一、数据前几行示例：

501
488
504
578
545
632
728
725
585

二、 前期准备

1.清理工作空间，设置工作目录并加载数据

rm(list=ls())
setwd("C:/Users/Lixue/Downloads") 
data <- scan("hotel.txt")

2. 数据准备与基础参数设置

T <- length(data)  # 数据的长度，即时间序列的总期数
s <- 12            # 数据的频率，每年有12个周期（即每月一个观察值）
years <- T/s       # 数据代表的年数
startyear <- 1     # 设置数据的起始年份
h <- 12            # 预测长度，设为12个月
fstart <- c(15,1)  # 预测开始的时间点：第15年的1月

3. 将数据转换为时间序列并绘制图形

data <- ts(data, start=c(startyear,1), frequency=s)
#将读取的数据转换为时间序列格式，起始时间为第1年的1月，频率为12（月度数据）
plot(data)

【ts()参数解释】

• start （）用于定义数据的起始时间。start=c(year, period) 是用来指定时间序列数据从哪一年、哪一周期开始的。start=c(1,1) 意思是时间序列数据从第1年的第1期（即第1月）开始。
• frequency=（s），R可以识别出这是一种月度时间序列。数据只有一列表示这是一个单变量时间序列（每个时间点上只有一个观测值，例如每个月的一个房间均价）。虽然数据是单列的，但R仍然需要知道数据的时间维度（比如起始年份和每年的周期数）。

4. 对数变换与趋势分析

y <- log(data)
#对数据进行对数变换，用于平滑季节性波动，并稳定数据的方差
plot(y)

三、构建趋势模型并拟合

1.模型构建与结果

t <- (1:T)  # 构建时间回归变量
model_l <- lm(y ~ t)  # 进行线性回归，以时间 t 作为自变量，y（对数变换后的数据）作为因变量，估计趋势线
summary(model_l)  # 查看回归模型的详细结果

Call:
lm(formula = y ~ t)
Residuals:
     Min       1Q   Median       3Q      Max 
-0.18664 -0.11001 -0.02185  0.04699  0.30596 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.3299787  0.0212156  298.36   <2e-16 ***
t           0.0027681  0.0002178   12.71   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1369 on 166 degrees of freedom
Multiple R-squared:  0.4933,	Adjusted R-squared:  0.4902 
F-statistic: 161.6 on 1 and 166 DF,  p-value: < 2.2e-16

2. 绘制趋势线

l_trend <- y - model_l$residuals
#获取回归模型的残差，并从原数据中减去残差，得到拟合的趋势线
plot(y, main="Add title here", ylab="Add title to y axis here")
lines(l_trend, col="green")
#将趋势线添加到原始数据的图形上，使用绿色显示

3. 计算增长率

growth <- exp(model_l$coefficients[2])  # 计算线性趋势的增长率
#coefficients[2]：获取回归系数中与时间 t 对应的斜率，表示每期的对数增长率
#exp()：计算实际的增长率（指数变换）
growth <- (growth - 1) * s * 100  # 将增长率转换为年化增长率

> growth
      t 
3.32638

四、去除趋势并分析季节性

• 【detrended <- y - l_trend 】

目的是去除时间序列中的趋势成分，得到去趋势后的时间序列，以便更清楚地观察时间序列中的季节性或其他周期性模式。

具体来说：

• 趋势反映的是时间序列数据中的长期变化，比如逐渐增加或减少的趋势。
• 季节性反映的是时间序列中每年或每月固定的周期性变化。

通过去除趋势，剩下的部分就是时间序列中的季节性和随机波动，去趋势后的数据可以更好地观测季节性模式或周期性变化

detrended <- y - l_trend
plot(detrended, ylab="add label here", main="Add Title Here")
abline(h=0, col="red")
#去除趋势成分，得到季节性波动的数据，并绘制图形，观察数据中纯粹的季节性因素

2. 计算季节性因素

• 【SeasonalFactor <- apply(seasonal.matrix, 2, mean)】

计算各个月份的平均季节性因素。

•  seasonal.matrix：矩阵每一列代表特定月份的季节性因素，每一行代表不同年份同一月份的观测值。比如，第一列可能代表每年的1月，第二列代表每年的2月，以此类推。可以通过时间序列的季节性分解得到。
• apply() 函数：主要用于在矩阵的行或列上进行操作。`apply(seasonal.matrix, 2, mean)` 表示对 `seasonal.matrix` 的每一列（参数 `2` 表示沿着列方向操作）应用 `mean()` 函数（每一列的平均值，也就是每个月在多年的数据中具有的平均季节性因素）。

seasonal.matrix <- matrix(detrended, ncol=s, byrow=TRUE)
#将去趋势后的数据按照月份（12个月一组）重组为矩阵
SeasonalFactor <- apply(seasonal.matrix, 2, mean) 
#按列计算每个月的平均值，得到季节性因素
SeasonalFactor <- ts(SeasonalFactor, frequency=s)
plot(SeasonalFactor)
abline(h=0, col="red")

3. 创建季节性哑变量

library("forecast") #forecast包，用于时间序列建模、预测和可视化
M <- seasonaldummy(y) #创建月度季节性哑变量
modelTS <- tslm(y ~ trend + M) # 创建线性模型
summary(modelTS)

• 【M <- seasonaldummy(y):】

创建月度季节性哑变量。seasonaldummy() 函数生成一组季节性哑变量（M1 到 M11），每个变量代表一个月份，用于捕捉每月的季节性变化。

M前十行示例：
       Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov
  [1,]   1   0   0   0   0   0   0   0   0   0   0
  [2,]   0   1   0   0   0   0   0   0   0   0   0
  [3,]   0   0   1   0   0   0   0   0   0   0   0
  [4,]   0   0   0   1   0   0   0   0   0   0   0
  [5,]   0   0   0   0   1   0   0   0   0   0   0
  [6,]   0   0   0   0   0   1   0   0   0   0   0
  [7,]   0   0   0   0   0   0   1   0   0   0   0
  [8,]   0   0   0   0   0   0   0   1   0   0   0
  [9,]   0   0   0   0   0   0   0   0   1   0   0
 [10,]   0   0   0   0   0   0   0   0   0   1   0

模型结果：

Call:
tslm(formula = y ~ trend + M)
Residuals:
      Min        1Q    Median        3Q       Max 
-0.059088 -0.013226  0.001383  0.014155  0.049916 
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.288e+00  6.427e-03 978.262  < 2e-16 ***
trend        2.725e-03  3.379e-05  80.653  < 2e-16 ***
MJan        -4.161e-02  8.016e-03  -5.190 6.50e-07 ***
MFeb        -1.121e-01  8.015e-03 -13.984  < 2e-16 ***
MMar        -8.446e-02  8.013e-03 -10.540  < 2e-16 ***
MApr         3.983e-02  8.012e-03   4.972 1.74e-06 ***
MMay         2.040e-02  8.011e-03   2.546   0.0119 *  
MJun         1.469e-01  8.010e-03  18.340  < 2e-16 ***
MJul         2.890e-01  8.009e-03  36.085  < 2e-16 ***
MAug         3.112e-01  8.009e-03  38.857  < 2e-16 ***
MSep         5.599e-02  8.008e-03   6.991 7.63e-11 ***
MOct         3.954e-02  8.008e-03   4.938 2.02e-06 ***
MNov        -1.122e-01  8.008e-03 -14.013  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.02119 on 155 degrees of freedom
Multiple R-squared:  0.9887,	Adjusted R-squared:  0.9878 
F-statistic:  1127 on 12 and 155 DF,  p-value: < 2.2e-16

结果显示，模型显著 (F=1127，p-value<0.0001)，线性趋势t和每个月度季节性哑变量均显著 (每个t-statistic的p-value小于0.05) 

4. 预测拟合值与可视化

fitTS <- y - modelTS$residuals
#从时间序列 y 中减去模型的残差，得到拟合的趋势和季节性部分 fitTS
plot(y, main="Add title here", ylab="Add title to y axis here")
lines(fitTS, col="green")
#绘制原始时间序列 y 并叠加模型拟合值 fitTS（绿色线）

五、两种方法进行预测

1. 使用 `season` 简化模型构建

如果数据具有线性趋势，一个简单的模型估计方式为：

• tslm(y ~ trend + season): 使用 season 来自动处理季节性，season 自动生成季节性因素，而不需要手动创建哑变量
• 这种方法自动将12个月中的1月排除（默认基线），而不是12月

model.alternative <- tslm(y ~ trend + season)
summary(model.alternative)
SDummyModelForecast <- forecast(model.alternative, h=h)
#使用模型进行 h 步预测
plot(SDummyModelForecast)

2. 手动预测下一年的12个月

对未来时间段进行预测，并将预测值从对数空间（log-transformed）转换回原始尺度。

1. 【tf <- ((T+1):(T+h))】

• 定义时间变量 tf，用于表示预测期间的时间点。
• T: 表示当前时间序列数据的长度，即观测数据的最后一个时间点。假设 T 是最后一个时间点，h 是预测的期数（比如预测未来12个月，则 h=12）。
• T+1:(T+h): 生成一个从 T+1 到 T+h 的序列，表示未来 h 个时间点。

2. 【MF <- I(seasonaldummy(y, h))】

• seasonaldummy(y, h): 生成未来 h 个时间点的季节性哑变量。
  • seasonaldummy(): forecast 包中的一个函数，根据时间序列 y 的频率（如月度或季度数据）生成季节性哑变量。
  • h: 表示生成 h 个未来时间点的季节性哑变量。假设 y 是月度数据，并且有12个月的周期（频率），那么 seasonaldummy(y, h) 将为未来12个月生成12组季节性哑变量。
• I(): 用于保护表达式防止被解释或进一步处理。这里的作用是确保生成的季节性哑变量 MF 被视为一个单独的对象，而不是进行公式解释。

结果：MF 包含未来 h 个时间点的季节性哑变量，形式类似于已经训练的模型中的月份哑变量。

tf <- ((T+1):(T+h)) #生成预测期的时间点
MF <- I(seasonaldummy(y, h)) 
Z <- data.frame(trend=tf, M=MF)
#创建数据框 Z ，包含了模型 modelTS 所需的所有自变量，包括趋势（时间序列的时间点）和未来时间点的季节性哑变量
prediction <- predict(modelTS, Z, interval="predict") #指定要计算预测区间，即置信区间。这个选项会返回三个值：点预测（即预测值）和上下两个置信区间（通常是95%的置信区间）
prediction <- exp(prediction)
#因为原始数据是对数变换的，预测值需要通过指数变换还原到原始尺度
prediction <- ts(prediction, frequency=s, start=fstart)
#prediction 被转换为一个时间序列对象

> prediction
             fit       lwr       upr
Jan 15  817.7227  782.7307  854.2790
Feb 15  764.1587  731.4588  798.3205
Mar 15  787.7030  753.9956  822.9173
Apr 15  894.3870  856.1143  934.3706
May 15  879.5636  841.9253  918.8845
Jun 15 1000.9106  958.0796 1045.6563
Jul 15 1156.9054 1107.3990 1208.6248
Aug 15 1186.0707 1135.3163 1239.0940
Sep 15  921.4229  881.9934  962.6151
Oct 15  908.8690  869.9767  949.5000
Nov 15  783.0270  749.5197  818.0322
Dec 15  878.4050  840.8162  917.6741

预测结果显示，在95%置信区间内，月度平均宾馆房间占有量在期间169（第15年的一月）处于783.73-854.28之间。

3. 生成预测图（含置信区间）

• xlim=c(startyear, 16)：指定X轴的取值范围。向量，表示X轴（时间轴）的取值范围。c(startyear, 16) 表示从 startyear 到16年（可以是年份，或者时间序列的相对时间点）
• ylim=c(min(data), 1.2*max(data))：指定Y轴的取值范围。

plot(data, xlim=c(startyear, 16), ylim=c(min(data), 1.2*max(data)), main="Add title here", ylab="Add title to y axis here")
fitTS <- exp(fitTS)
lines(fitTS, col="green")
lines(prediction[,1], col="blue") #点预测值
lines(prediction[,2], col="red") #下置信区间
lines(prediction[,3], col="red") #上置信区间

4. 残差分析

• 【dwtest()】

进行Durbin-Watson检验，检查残差中的自相关性。

Durbin-Watson统计量的范围：

• Durbin-Watson统计量 (DW) 的值介于 0 和 4 之间：
  • DW≈2：没有自相关，残差是独立的。
  • DW<2：正自相关，表明相邻观测值之间的残差倾向于同方向变化。
  • DW>2：负自相关，表明相邻观测值之间的残差倾向于相反方向变化。

plot(modelTS$residuals)
dwtest(modelTS)

Durbin-Watson test
data:  modelTS
DW = 1.19, p-value = 1.707e-07
alternative hypothesis: true autocorrelation is greater than 0

结果显示，模型可能存在自相关问题（Durbin-Watson统计量 = 1.19），需要进一步改进。

六、使用三角函数模型建模季节性

1.参数设置

• sin() 和 cos(): 使用三角函数捕捉季节性周期。
  • st1 和 st2 表示一年中的季节性波动（基本周期）。
  • st3 和 st4 表示较高频的季节性波动（两倍频率）。
• lm(): 构建线性模型，将 y 拟合到趋势 t 和三角函数 st1 到 st4 上。

t <- (1:T)
st1 <- sin(2 * pi * t / s)
st2 <- cos(2 * pi * t / s)
st3 <- sin(4 * pi * t / s)
st4 <- cos(4 * pi * t / s)
model.tri <- lm(y ~ t + st1 + st2 + st3 + st4)
summary(model.tri)

Call:
lm(formula = y ~ t + st1 + st2 + st3 + st4)
Residuals:
      Min        1Q    Median        3Q       Max 
-0.100881 -0.045403 -0.002315  0.047136  0.122679 
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.3327488  0.0086962 728.222  < 2e-16 ***
t            0.0027354  0.0000893  30.630  < 2e-16 ***
st1         -0.1009242  0.0061215 -16.487  < 2e-16 ***
st2         -0.1266548  0.0061131 -20.719  < 2e-16 ***
st3          0.0662648  0.0061144  10.838  < 2e-16 ***
st4          0.0189788  0.0061131   3.105  0.00225 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.05602 on 162 degrees of freedom
Multiple R-squared:  0.9172,	Adjusted R-squared:  0.9146 
F-statistic: 358.7 on 5 and 162 DF,  p-value: < 2.2e-16

结果显示，模型显著(F= 358.7 and  p-value < 0.001) ，线性趋势和三角函数表示的季节性波动也均显著 (每个t-statistic的p-value均小于0.05). R方0.9172 ，标准误0.05602。

2.预测并可视化三角函数模型的结果

使用三角函数模型进行预测，生成未来 h 期的预测区间，并通过指数变换恢复预测值。

triZ <- data.frame(t=tf, st1=sin(2*pi*tf/s), st2=cos(2*pi*tf/s), st3=sin(4*pi*tf/s), st4=cos(4*pi*tf/s))
triprediction <- predict(model.tri, triZ, interval="predict")
triprediction <- exp(triprediction)
triprediction <- ts(triprediction, frequency=s, start=c(years+1, 1))

> triprediction
             fit       lwr       upr
Jan 15  813.8318  726.6624  911.4579
Feb 15  808.2960  721.6543  905.3399
Mar 15  796.7918  711.3655  892.4768
Apr 15  822.4436  734.2854  921.1861
May 15  913.5094  815.6099 1023.1601
Jun 15 1047.6827  935.3979 1173.4462
Jul 15 1139.6342 1017.4693 1276.4670
Aug 15 1110.7615  991.6851 1244.1360
Sep 15  991.1383  884.9077 1110.1215
Oct 15  877.2949  783.2853  982.5875
Nov 15  824.8963  736.4836  923.9226
Dec 15  826.6974  738.0225  926.0268

3.比较预测结果

比较三角函数模型的预测结果与哑变量模型的预测结果，通过观察预测区间长度，可以看出三角函数模型预测区间更宽。

plot(data, xlim=c(startyear, 16), ylim=c(min(data), 1.2*max(data)), main="Add title here", ylab="Add title to y axis here")
lines(modelfit.tri, col="green")
lines(triprediction[,1], col="blue")
lines(triprediction[,2], col="red")
lines(triprediction[,3], col="red")

4.残差分析

plot(model.tri$residuals)
dwtest(model.tri)

Durbin-Watson test
data:  model.tri
DW = 2.6303, p-value = 0.9999
alternative hypothesis: true autocorrelation is greater than 0

• DW=2.6303：统计量接近3，表明残差有负自相关的趋势（即相邻的残差更可能朝相反的方向波动）。
• p-value = 0.9999：检验的p值，用来评估零假设（即没有自相关）是否成立。极高的p值（接近1）表明我们不能拒绝零假设，即模型残差没有显著的自相关性。说明模型在解释时间序列数据的方面表现良好，残差是随机的且独立分布的。