R语言时间滞后模型 r语言时间序列建模_拟合

1. 时间序列概述

按照时间的顺序把随机事件变化发展的过程记录下来就构成了一个时间序列

对时间序列进行观察、研究，找寻它变化发展的规律，预测它将来的走势就是时间序列分析

时间序列建模基本步骤：

R语言时间滞后模型 r语言时间序列建模_r语言_02

解释建模的基本步骤：

  • 通过read.table()收集数据，ts()绘制时序图
  • 根据观察时序图以及白噪声检验Box.test()，进行平稳性判别的检验
  • 若得到平稳的非白噪声序列，则进行模式识别
  • 画自相关图和非自相关图，根据两图的结尾性和拖尾性进行AR、MA、ARMA的模式识别
  • 对识别后模式中的位置参数进行参数估计arima()
  • 模型检验分为：①残差的白噪声检验；②过度拟合检验pt()
  • 模型检验通过则进行模型优化，否则重新进行模式识别
  • 模型优化中得到AIC和BIC值，进行模型的优化
  • 然后进行预测与控制

2. 平稳性判别

2.1 时序图绘制

时序图是一张二维平面图，一般横坐标表示时间，纵坐标表示序列取值，观察时序图，我们能获得序列值的趋势和走向

时序图的绘制用函数ts()

#将美国爱荷华州杜比克 (Dubic) 市 144 个月的平均气温 (单位：
#°F) 按时间顺序记录下来, 就得到长度为 144 的观察值序列
 > t <- scan("E:/DATA/CHAP1/data1.2.txt")
 > t <- ts(t, start=c(1964,1),frequency=12)
 > plot(t,type="o",xlab="年份",ylab="气温",col=4)

R语言时间滞后模型 r语言时间序列建模_拟合_03

根据平稳性的定义，平稳时间序列的均值和方差均为常数，因此平稳时间序列的时序图应该围绕一条水平线上下波动，而且波动的范围有界

如果序列时序图显示出了明显的趋势性或周期性，那么它通常不是平稳的时间序列

根据这个性质，许多时间序列通过时序图就可看出它的非平稳性

2.2 数据的白噪声检验

当识别出一个序列是平稳时间序列之后，我们需要进一步分析该序列是否为纯随机序列，因为如果是纯随机序列的话，那么意味着序列值之间没有相关关系，该序列就成为所谓的无记忆序列,即过去的行为对将来的发展没有丝毫影响，这样的序列从统计分析的角度而言无任何研究的意义

纯随机性检验也叫白噪声检验

在R中使用函数 Box.test() 进行白噪声检验

Box.test(x, type=, lag= )
  • x：变量名, 可以是数值向量, 也可以是一元时间序列名
  • type：检验统计量类型
#计算例 1.20 中白噪声序列分别延迟 6 期和 12 期的 统计量的值, 并判断该序列的随机性 (α = 0.05)
> Box.test(white_noise,lag=6)
 Box-Pierce test
 data: white_noise
 X-squared = 2.0304, df = 6, p-value = 0.9169
> Box.test(white_noise,lag=12)
 Box-Pierce test
 data: white_noise
 X-squared = 6.838, df = 12, p-value = 0.8681

2.3 平稳性判别

得到平稳的、非白噪声的时间序列后进行模式的识别

3. 模型识别

3.1 自相关图的绘制

在R中绘制自相关图使用函数acf( )：这个函数的命令格式为：

acf(x, lag)
  • x：是时间序列数据构成的向量
  • lag：是延迟的阶数。若用户不特殊指定的话，系统会根据序列长度
# 绘制2014年7月至2017年5月北京市商品住宅施工面积累计值 的时序图, 并添加辅助线来比较
> w <- read.csv("E:/DATA/CHAP1/Beijing commodity 
+housing.csv", header=T)
> z <- w$CCA
> y <- na.spline(z)
> CCA <- ts(y,start=c(2014,7),frequency = 12)
> plot(CCA,type="o",col=1)
> abline(v=2014.9,lty=2, col=1,lwd=2)
> abline(v=c(2015.1,2016.1),lty=2,col=1)
> abline(h=c(6261.21,5857.61), lty=3, col=1,lwd=2)
# 绘制 2014 年 7 月至 2017 年 5 月北京市商品住宅施工面积累计值的自相关图
> acf(CCA) 
#虚线为自相关函数的 2 倍标准差位置. 一般地, 如果悬垂线夹在两条虚线之间,那么可认为此时自相关函数非常接近于零

R语言时间滞后模型 r语言时间序列建模_时间序列_04

3.2 偏自相关的绘制

偏自相关图用pacf()来绘制

R语言时间滞后模型 r语言时间序列建模_拟合_05

3.3 模式的识别

观察自相关图和偏自相关图的图像的截尾性和拖尾性

R语言时间滞后模型 r语言时间序列建模_r语言_06

4. 序列的拟合

4.1 filter()拟合函数

filter( ) 函数可以直接拟合 AR 序列 (不论是否平稳) 和移动平均序列 (MA序列)。filter( ) 函数的命令格式为：

filter(a, filter= , method= , circular= )
  • a：随机波动序列的变量名
  • filter：指定模型系数
  • method：指定拟合的是AR模型还是MA模型：
  • method=“recursive” 为 AR 模型
  • method=“convolution” 为 MA 模型
  • circular：拟合 MA 模型时专用的一个选项, circular=T 可以避免NA数据出现

4.1 arima.sim拟合函数

arima.sim( ) 函数可以拟合平稳 AR 序列以及MA 序列、平稳自回归移动平均模型（ARMA）序列和 ARIMA序列

arima.sim(n, list(ar= ,ma= ,order= ),sd= )
  • n: 拟合序列的长度
  • list: 指定具体模型参数
  • sd: 指定序列的标准差, 如不特殊指定, 系统默认 sd=1

R语言时间滞后模型 r语言时间序列建模_时序图_07

R语言时间滞后模型 r语言时间序列建模_时序图_08

5. 参数估计

模式识别完成后，要对模型中未知参数进行估计

在R语言中, 参数估计可通过调用函数arima完成

arima(x, order=, include.mean=, method= )
  • x：要进行模型拟合的序列名
  • order：指定模型阶数。order=c(p,d,q)，其中p为自回归阶数; q 为移动平均阶数；d 为差分阶数，差分阶数后面章节才会用到，在本章取d = 0
  • include.mean：决定是否包含常数项. 如果 include.mean=T, 那么需要拟合常数项, 这是系统默认设置; 如果 include.mean=F, 那么意味着不需要拟合常数项
  • method：指定参数估计方法。如果 method=“CSS-ML”，那么指定参数估计方法是条件最小二乘和极大似然估计混合方法，这是系统默认设置；如果method=“ML”，那么指定参数估计方法是极大似然估计法；如果method=“CSS”，那么指定参数估计方法是条件最小二乘估计法
#确定 2016 年 1 月至 2017 年 6 月青海省居民消费价格指数序列拟合模型的口径(即对该序列的未知参数进行估计)
> x <-
read.table("E:/DATA/CHAP3/cpi.csv",sep=",",header=T)
> QHCPI <- ts(x$QHCPI,start=1)
> QHCPI.fix <- 
arima(QHCPI,order=c(2,0,0),method="ML")
> QHCPI.fix
Call:
arima(x = QHCPI, order = c(2, 0, 0), method = "ML")
Coefficients:
       ar1       ar2     intercept
      1.0233   -0.5067    101.5002
s.e.  0.1975    0.1928     0.2245
sigma^2 estimated as 0.2123: log likelihood = -12.2, aic = 32.4

6. 模型检验

6.1 残差检验

如果模型被正确识别，参数估计足够精确，那么残差应该具有白噪声的性质，即残差序列应表现出独立、同分布、零均值和相同标准差的性质。

反之，如果残差序列为非白噪声序列，那就意味着残差序列中还残留着相关信息未被提取，说明拟合模型不够有效，需要重新选择其他模型进行拟合。

因此, 残差的检验指的就是残差序列的白噪声检验

白噪声检验参考上述2.2

6.2 过度拟合检验

在模型诊断中，需要特别引起注意的问题是，由过度拟合而产生的参数冗余问题

pt(t, df=, lower.tail= )
  • t：t统计量的值
  • df：自由度
  • lower.tail：确定计算概论的放心
#对 2016 年 1 月至 2017 年 6 月青海省居民消费价格指数序列拟合模型的参数进行显著性检验
> x <- read.table("E:/DATA/CHAP3/cpi.csv",sep=",",header=T)
> QHCPI <- ts(x$QHCPI,start=1)
> QHCPI.fix <- arima(QHCPI,order=c(2,0,0),method="ML")
> QHCPI.fix
Call:
 arima(x = QHCPI, order = c(2, 0, 0), method = "ML")
 Coefficients:
 ar1 ar2 intercept
 1.0233 -0.5067 101.5002
s.e. 0.1975 0.1928 0.2245
sigma^2 estimated as 0.2123: log likelihood = -12.2, aic = 32.4

> # ar1 系数的显著性检验
> t1 <- 1.0233/0.1975
> pt(t1,df=15,lower.tail=F)
[1] 5.584387e-05
> # ar2 系数的显著性检验
> t2 <- -0.5067/0.1928 
> pt(t2,df=15,lower.tail=T)
[1] 0.009501987
> # 常数的显著性检验
> t0 <- 101.5002/0.2245 
> pt(t0,df=15,lower.tail=F)
[1] 9.938772e-33
# 三个系数检验的 p 值均小于 0.05, 故三个系数均显著非零

7. 模型优化

由于样本的随机性和定阶过程很大程度上依赖于分析人员的主观判断，所以在模型识别时，可能就会有若干个备选模型符合条件，而且有时会出现多个模型都通过了检验，那么现在的问题是，到底哪个模型最有效呢

一个好的拟合模型应该是在兼顾考虑拟合精度和未知参数的个数下，从中选择最优的配置，基于此AIC函数被提出

AIC准则认为，使得AIC函数达到最小的模型是最优模型

AIC准则为模型选择提供了重要依据，但是AIC准则也有不足之处

为了弥补AIC准则的不足，1976年提出了BIC准则

BIC准则规定，使得BIC函数达到最小的模型是最优模型

总之：AIC和BIC越小，则模型越优

在R语言中，通过调用程序包forecast中的函数Arima( )也可以进行参数估计，其命令格式和参数使用方式与函数arima( ) 的一样。所不同的是, 函数Arima( )的运行结果会同时给出AIC、BIC的值

> QHCPI.fix2$aic
 [1] 32.79913
> QHCPI.fix2$bic
 [1] 37.25099

8. 序列预测

所谓预测 (forecast)，就是根据现在与过去的随机序列的样本取值，对未来某个时刻序列值进行估计

在R中，我们可以通过调用函数forecast()来进行预测

forecast(object, h=, level= )
  • object：拟合信息文件名
  • h：预测期数
  • level：置信区间的置信水平