R语言时间序列学习+代码记录2：Modeling cycles - AR( ), MA( ), and ARMA( )

数据来源主要为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

一、数据示例：

a quarterly seasonally adjusted index of Canadian employment, 1962.1 - 1993.4
num [1:136]
83.09025 
82.79963 
84.63444 
85.37746 
86.19760 
86.57884 
88.04972
87.92493 
88.46513 
88.39846
... ...

二、 前期准备

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

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

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

T <- length(data)  # 数据的长度，即时间序列的总期数
s <- 4            # 数据的频率，季度数据
years <- T/s       # 数据代表的年数
styear <- 1962     # 设置数据的起始年份（Q1 1962）
h <- 12            # 预测长度
fstart <- c(1994,1)  # 预测开始的时间点（Q1 1994）
years <- T/s       #计算数据覆盖年数（每年均有Q1~Q4）

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

data <- ts(data, frequency=s, start=datastart) 
#将读取的数据转换为时间序列格式，起始时间为第1962年的Q1，频率为4（季度数据）
plot(data)

由图像：序列没有显示趋势和季节性 （数据已经过季节性调整），然而，数据中的就业情况具有强烈的周期性特征。就业情况似乎与时间序列高度相关——在经济繁荣时期就业率较高，而在经济衰退时就业率较低。

 4. ACF, PACF函数探索自相关性

y <- data
acf(y)
pacf(y)

• acf(y)：绘制自相关函数 (ACF)，显示数据的滞后自相关情况。
• pacf(y)：绘制偏自相关函数 (PACF)，用于观察各滞后阶的自相关性。

ACF 和 PACF 是用于帮助确定 ARMA 模型阶数的关键工具。MA 模型的阶数通常可以从 ACF 图中看到，ACF 会在某个滞后期截断（即滞后q），而 AR 模型的阶数则可以从 PACF 中截断的滞后期确定。

ACF:

PACF：

三、拟合移动平均模型 (MA(q))

1.自相关函数 (ACF)：

• 滞后 (Lag)：某个时刻的数据与之前某个时间点的数据之间的差距。例如，滞后1表示当前时间点的数据与前一个时间点的数据之间的关系。
• 自相关系数 (Autocorrelation Coefficient)：测量当前时间点的数据与滞后期数据之间的线性相关性，范围通常在 [-1, 1] 之间。
• 显著性界限 (Significance Boundaries)：在 ACF 图中，通常会绘制水平线，表示零相关性 ± 2 个标准误的范围。如果自相关系数超出了这个范围，说明滞后期的自相关性显著。

在移动平均模型 (MA) 中，自相关函数 (ACF) 的重要作用在于识别模型的阶数。对于一个 MA(q) 模型，ACF 会在滞后期 q 之后 截断 (cut-off)，即在 q 之后的滞后期自相关系数迅速衰减至零。因此，通过观察 ACF 中截断的位置，可以帮助确定 MA 模型的阶数 q（MA（q））。

2. MA模型

R语言中，可以使用arima函数拟合MA(q)模型，参数order应设置为c(0,0,q)。由图像，数据中存在显著的自相关性，ACF显示有11个显著的峰值，因此使用q=11的MA模型

• arima() 函数：拟合 ARIMA 模型。其中，order=c(0,0,11) 表示拟合一个 MA(11) 模型：
  • 第一个参数 0 代表 AR(自回归) 部分的阶数为0。
  • 第二个参数 0 代表差分次数，即假设数据是平稳的（不需要差分）。
  • 第三个参数 11 代表 MA(移动平均) 部分的阶数为11。

model.ma <- arima(y,order=c(0,0,11))
model.ma

Call:
arima(x = y, order = c(0, 0, 11))

Coefficients:
         ma1     ma2     ma3     ma4     ma5     ma6     ma7     ma8
      1.4758  1.7407  1.8294  1.8091  1.6638  1.4703  1.1177  0.8466
s.e.  0.0848  0.1570  0.2286  0.2798  0.2965  0.2762  0.2319  0.1972
         ma9    ma10    ma11  intercept
      0.7216  0.5742  0.2750    99.5410
s.e.  0.1867  0.1604  0.0882     1.7405

sigma^2 estimated as 2.08:  log likelihood = -245.15,  aic = 516.29

 由于模型阶数较高，在训练时丢失了前 11 个数据点，这是高阶模型的一个缺点。出于模型复杂性和数据可用性的考虑，与参考书 (Diebold 的《Elements of Forecasting》) 保持一致，尝试一个较低阶的模型 MA(4)。

• MA(4)在拟合过程中只使用前 4 个滞后期来解释时间序列中的依赖关系，从而保留了更多的原始数据点，避免了 MA(11) 模型中可能出现的过度拟合问题。

model.ma <- arima(y, order=c(0,0,4))
model.ma

Call:
arima(x = y, order = c(0, 0, 4))

Coefficients:
         ma1     ma2     ma3     ma4  intercept
      1.6724  1.7995  1.3046  0.5492    99.9243
s.e.  0.0712  0.1150  0.1145  0.0649     1.0278

sigma^2 estimated as 3.685:  log likelihood = -283.63,  aic = 579.26

3. 检查残差

• plot(model.ma$residuals) ：绘制MA(q) 模型的残差图。显示模型所不能解释的部分信息，应该接近白噪声（即无序且无自相关性）。
• acf(model.ma$residuals) ：计算并绘制残差的自相关函数图。

如果 MA(q) 模型能够有效地捕捉数据的动态特征，残差应该呈现出白噪声特征：残差的 ACF 图应该显示出所有滞后期的自相关系数都接近于零，且在显著性界限之内。

plot(model.ma$residuals)
acf(model.ma$residuals)

MA（11） ：

MA（4）： 

由图：MA(4) 模型的残差 ACF 图仍然存在显著的自相关性，表明模型未能完全捕捉到时间序列数据中的自相关特征，即模型未能充分拟合数据，遗漏了一些周期性或趋势性。可能需要更高阶的模型或者考虑使用其他模型类型（例如 AR 模型）。

四、拟合自回归模型 (AR（p）)

1. 偏自相关函数 (PACF)：

• 直接相关性 (Direct Correlation)：排除了中间滞后期的影响后，滞后期 k 与当前时间点之间的直接相关性。
• 截断性 (Cut-off)：与 ACF 不同，PACF 在 AR(p) 中，在滞后期 p 之后截断（即在 p 之后的滞后期偏自相关系数迅速衰减至零）。

PACF 衡量的是在控制了所有中间滞后期的影响后，不同滞后期之间的直接相关性。具体来说，PACF 是当前时间点与滞后期 k 之间的相关性，排除了所有比 k 小的滞后期对两者相关性的影响。

2. 自回归（AR）模型

描述一个变量当前值与其过去值之间的线性关系。当前时间点的值由之前的自身数据（滞后数据）线性回归而得出。 AR(p) 模型中 p 代表模型中所使用的滞后期数。

PACF 图中，滞后期 2 后出现截断，尝试拟合 AR(2) 模型，即假设当前时间点的数据与前两个滞后期有较强的相关性。

•  arima() 函数：拟合 AR 模型。order=c(2,0,0) 表示拟合一个 AR(2) 模型：
  • 第一个参数 2 代表自回归部分的阶数 p=2。
  • 第二个参数 0 代表差分次数，即假设数据是平稳的（不需要差分）。
  • 第三个参数 0 代表移动平均部分的阶数为 0，即不包含 MA 组件。

model.ar <- arima(y, order=c(2,0,0))
model.ar

Call:
arima(x = y, order = c(2, 0, 0))

Coefficients:
         ar1      ar2  intercept
      1.4505  -0.4763    97.4982
s.e.  0.0749   0.0762     4.3936

sigma^2 estimated as 2.022:  log likelihood = -242.79,  aic = 493.57

3. 检查残差

• plot(model.ar$residuals) ：绘制 AR(2) 模型的残差图，检查模型是否很好地拟合了数据。如果残差看起来像白噪声（即没有明显的模式或结构），则表明模型适用。
• acf(model.ar$residuals) ：计算并绘制残差的自相关函数图，进一步验证残差的随机性。

如果模型适当地拟合了数据，则残差的 ACF 应该在所有滞后期上都接近零，表明残差没有自相关性。

由图：结果良好，残差表现为白噪声，这意味着 AR(2) 模型能够很好地解释数据中的自相关特性。

五、自回归移动平均模型（ARMA(p, q) ）

1.赤池信息准则 (AIC)

AIC 是用于选择模型的一个标准，计算公式为：AIC=2k−2ln(L)

• 模型的似然函数 L：表示模型的拟合好坏。
• 模型中的参数数量 k
• AIC 值越小，表示模型的拟合效果越好，同时考虑了模型的复杂度。

（1）使用 auto.arima() 函数自动选择最佳模型

函数根据给定的数据自动搜索最佳的 ARIMA 或 ARMA 模型。因为指定了 d=0 和 seasonal = FALSE，它会在 ARMA 模型中搜索最佳的 ARMA(p, q) 组合。

参数设置：

• d=0：指定差分次数为 0，假设数据已经是平稳的（不需要差分）。
• seasonal = FALSE：表示不考虑季节性成分。
• ic = "aic"：指定使用赤池信息准则 (AIC) 作为模型选择的标准。AIC 值越小，模型越好。
• stepwise = FALSE：关闭逐步选择（即试探更多的模型组合）。
• approximation = FALSE：不使用近似估计，确保估计精确。
• trace = TRUE：输出所有尝试的模型和对应的 AIC 值，方便选择最优模型。

require(forecast)
model.arma <- auto.arima(y, d=0, seasonal = FALSE, ic = "aic", 
                         stepwise = FALSE, approximation = FALSE, 
                         trace = TRUE)

 ARIMA(0,0,0)           with zero mean     : 1642.012
 ARIMA(0,0,0)           with non-zero mean : 954.4594
 ARIMA(0,0,1)           with zero mean     : Inf
 ARIMA(0,0,1)           with non-zero mean : 795.8631
 ARIMA(0,0,2)           with zero mean     : Inf
 ARIMA(0,0,2)           with non-zero mean : 686.9732
 ARIMA(0,0,3)           with zero mean     : 1153.926
 ARIMA(0,0,3)           with non-zero mean : 624.4021
 ARIMA(0,0,4)           with zero mean     : Inf
 ARIMA(0,0,4)           with non-zero mean : 579.2624
 ARIMA(0,0,5)           with zero mean     : Inf
 ARIMA(0,0,5)           with non-zero mean : 562.6548
 ARIMA(1,0,0) with zero mean     : Inf
 ARIMA(1,0,0)           with non-zero mean : 525.9286
 ARIMA(1,0,1) with zero mean     : Inf
 ARIMA(1,0,1)           with non-zero mean : 502.0656
 ARIMA(1,0,2) with zero mean     : Inf
 ARIMA(1,0,2)           with non-zero mean : 497.5256
 ARIMA(1,0,3) with zero mean     : Inf
 ARIMA(1,0,3)           with non-zero mean : 498.4639
 ARIMA(1,0,4) with zero mean     : Inf
 ARIMA(1,0,4)           with non-zero mean : 499.1958
 ARIMA(2,0,0) with zero mean     : Inf
 ARIMA(2,0,0)           with non-zero mean : 493.5718
 ARIMA(2,0,1) with zero mean     : Inf
 ARIMA(2,0,1)           with non-zero mean : 494.8726
 ARIMA(2,0,2) with zero mean     : Inf
 ARIMA(2,0,2)           with non-zero mean : 496.6441
 ARIMA(2,0,3) with zero mean     : Inf
 ARIMA(2,0,3)           with non-zero mean : 497.8419
 ARIMA(3,0,0) with zero mean     : Inf
 ARIMA(3,0,0)           with non-zero mean : 494.9481
 ARIMA(3,0,1) with zero mean     : Inf
 ARIMA(3,0,1)           with non-zero mean : Inf
 ARIMA(3,0,2) with zero mean     : Inf
 ARIMA(3,0,2)           with non-zero mean : 498.7393
 ARIMA(4,0,0) with zero mean     : Inf
 ARIMA(4,0,0)           with non-zero mean : 496.8779
 ARIMA(4,0,1) with zero mean     : Inf
 ARIMA(4,0,1)           with non-zero mean : 498.8077
 ARIMA(5,0,0)           with zero mean     : Inf
 ARIMA(5,0,0)           with non-zero mean : 498.4932

 Best model: ARIMA(2,0,0)           with non-zero mean 

输出结果表明：

• 最佳选择：ARMA(2,0) 模型（即 AR(2) 模型）是 AIC 最小的模型。
• 次优模型：ARMA(2,1) 模型，尽管稍微复杂一些（包含一个 MA 部分），但其拟合效果并没有显著优于 AR(2) 模型。

2. ARMA 模型

ARMA 模型是将自回归 (AR) 模型和移动平均 (MA) 模型结合起来的一种时间序列模型，能够更灵活地捕捉数据中的动态特征。ARMA(p, q) 模型包含两个部分：

• AR(p)：代表自回归部分，用过去的 p 个时间点的值来解释当前值。
• MA(q)：代表移动平均部分，用当前及过去 q 个误差项来解释当前值。

（1）手动拟合 ARMA(2,1) 模型

使用 arima() 函数手动拟合 ARMA(2,1) 模型，其中 order=c(2,0,1) 表示：

• 2 是 AR 部分的阶数，表示使用前两个滞后值。
• 0 是差分次数，表示数据不做差分。
• 1 是 MA 部分的阶数，表示使用一个滞后误差项。

model.arma <- arima(y, order=c(2,0,1))
model.arma

Call:
arima(x = y, order = c(2, 0, 1))

Coefficients:
         ar1      ar2      ma1  intercept
      1.5745  -0.5987  -0.1612    97.9320
s.e.  0.1534   0.1522   0.1922     4.0145

sigma^2 estimated as 2.011:  log likelihood = -242.44,  aic = 494.87

结果表明：尽管 ARMA(2,1) 模型拟合效果不错，但与 AR(2) 模型相比并没有显著改进，因此 AR(2) 模型是更优选择。 

 3. 检查残差是否为白噪声

绘制残差的 自相关函数 (ACF) 和 偏自相关函数 (PACF) 图。理想情况下，残差应该呈现出白噪声特性，即没有显著的自相关性

总结：ACF 与 PACF

1.选择合适的模型阶数

根据模型类型不同，ACF 和 PACF 的图形特征各有不同：

• AR(p) 模型：

  • PACF 在滞后期 p 后截断。
  • ACF 逐渐衰减，可能以指数形式减小。

• MA(q) 模型：

  • ACF 在滞后期 q 后截断。
  • PACF 逐渐衰减。

• ARMA(p,q) 模型：

  • ACF 和 PACF 都逐渐衰减，没有明显的截断点。

通过比较 ACF 和 PACF 的图形特征，可以选择合适的模型：

• 如果 ACF 截断而 PACF 衰减，则适合 MA 模型。
• 如果 PACF 截断而 ACF 衰减，则适合 AR 模型。
• 如果 ACF 和 PACF 都逐渐衰减，则可以考虑 ARMA 模型。

2.残差图像：

（1）残差 ACF 图：

衡量不同滞后期的残差自相关性。理想情况下，若残差是白噪声，则所有滞后期的自相关系数应在蓝色虚线（显著性边界）内。

在滞后 k 处，若自相关系数显著高于显著性边界，则表明在该滞后期存在显著自相关性，说明模型的残差仍存在一些结构性关系，可能需要进一步改进，尤其是考虑自回归部分的修正。

（2）残差 PACF 图：

用于测量当前残差与不同滞后期残差之间的直接自相关性（即在控制其他滞后影响后的自相关性）。

若在滞后 k 处，偏自相关系数接近显著性界限，若自相关系数显著高于显著性边界，则表明可能还有某些小范围的自相关结构没有被模型捕捉。