Old John World

Welcome to my world

R语言-时间序列分析

发表于 更新于 分类于
本文字数： 6.3k 阅读时长 ≈ 12 分钟

本文介绍时间序列模型的基本r语言代码，包括平稳性和非平稳性的处理。

一个完整的Arma模型分析处理

arma模型平稳性判别

1. 时序图绘制

1
2
3
4
5
6
7
# 1. 时间序列的赋值
rain<-ts(data$price,start=2015)
# 时间序列的子序列
rain2<-window(rain,start=2016,end=2018)
# data表格的price列，开始是2015年
# 2. 绘制
plot()

3. 噪声检验（纯随机检验）

1
2
3
4
for(k in 1:2) print(Box.test(white_noise,lag=6*k,type="Ljung-Box"))#噪声检验
# type='Box-Pierce'是Q检验量
## 输出LB统计量和P值
## p值大于0.05则不能拒绝

结果分析：显示延迟6和12阶的LB统计量的P值小于显著性水平（α = 0.05)，所以可以判断该序列为平稳非白噪声序列。

DF/ADF平稳性检验

1
2
3
4
install.packages('aTSA')
library(aTSA)
adf.test(overshort,nlag=)
# P值为0.01小于0.05，所以该序列平稳

DF检验针对AR(1)，是ADF的特例。
结果分析：6中子类型，τ统计量的P值均显著大于显著性水平（α = 0.05)，因此可以判断，如果序列考虑如上6种结构之一提取确定性信息，则随机性部分都不能实现平稳。
所以序列是非平稳序列。
若存在小于显著性水平的，则是平稳序列。

2. 自相关图and偏自相关图（模型定阶）

1
2
3
4
# lag延迟阶数
# plot=True，只输出自相关图，不输出系数；False，反之。
acf(x,lag.max=,plot=)
pacf(x,lag.max=,plot=)

结果分析：
判断是截尾还是拖尾，然后分析模型属于哪一种。

自相关图是自相关系数，虚线是自相关系数两倍标准差的参考线。落在两倍标差外，则认为自相关系数很大，显著非零。

参数估计

1
2
3
4
x.fit<-arima(x,order=c(p,d,q),include.mean=,method=)
# order为模型阶数:自回归阶数、差分阶数、移动平均阶数。
# include.mean要不要包含均值
# method="CSS-ML"最小二乘和极大似然的混合 也可以单独使用一种

模型检验：残差分析

1
2
library(aTSA)
ts.diag(x.fit)

结果输出四图

ACF/PACF图：需要所有竖线全都在蓝色虚线内(95%CI)。

第三幅图：残差序列的白噪声检验，观察各阶延迟下的白噪声检验统计量的P值是否都显著的大于0.05。若都显著，则残差序列属于白噪音序列，拟合模型显著成立。
第四幅图：QQ图，正态分布假定的检验。如果点密集地分布在对角线左右，则该序列近似服从正态分布。

参数的显著性检验

检验每一个未知参数是否显著非零
目的：精简模型

1
2
3
4
5
6
7
# 1. 近似方法
x.fit
# 2. 精确方法
# 构造t统计量，调用pt函数求p值
t = abs(fit$coef)/sqrt(diag(fit2$var.coef))
pt(t,length(overshort)-length(fit1$coef),lower.tail=F)
# t统计量的值，自由度，参数估计值为正

模型优化

1
2
BIC(x)
AIC(x)

识别最优阶数

1
2
3
4
install.package('forecast')
library(forecast)
auto.arima(x,ic='bic')
auto.arima(x,ic='aic')# 默认

预测

1
2
3
library(forecast)
fore<-forecast::forecast(fit,h=10)#预测未来十个数据
plot(x.fore)

非平稳处理

线性拟合

1
2
3
4
5
6
7
8
9
10
x <-c()
# 40年
t <-c(1:40)
x.fit <- lm(x~t)
# 查看拟合信息
summary(x.fit)
# 画拟合图
x <- ts(x)
plot(x)
abline(lm(x~t),col=2)

曲线拟合

可用线性最小二乘

1
2
3
4
5
t1<-c(1:60)
## 转化
t2<-t1^2
x.fit1<-;m(x~t1+t2)
summary(x.fit1)

不可用线性最小二乘

1
2
3
4
5
6
7
x.fit2<-nls(x~a+b*t1+c*t1^2,start=list(a=1,b=1,c=1))
summary(x.fit2)
## 画图
y<-predict(x.fit2)#把nls函数得到的拟合值赋值给y
y<-ts(y,start=1949)
plot(x,type="p")
lines(y,col=2,lwd=2)

移动平均拟合

1
2
3
4
5
6
7
8
9
# 移动平均拟合
library(TTR)
x<-
# 选择期数 进行拟合
x.ma<-SMA(x,n=5)
# 画图
plot(x,tyoe='o')
lines(x.ma,col=2,lwd=2)

指数平滑

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
x<-
x.fit<-HoltWinters(x,alpha=, beta=, gamma=F,seasonal=)
# 简单指数平滑
x.fit<-HoltWinters(x,beta=F,gamma=F)
# Holt两参数平滑
x.fit<-HoltWinters(gamma=F)
# Holt'三参数平滑
x.fit<-HoltWinters()
# seasonal：季节和趋势的关系
# seasonal='additive';seasonal='multiplicative'
# 默认加法，seasonal='mult'乘法
#画图
plot(x.fit)
x.fore<-forecast(x.fit,h=10)
plot(x.fore)

季节指数、趋势效应、随机效应分解

1
2
3
4
5
6
7
8
9
10
x <- ts()
# 乘法模型
x.fit <- decompose(x,type='mult')
# 加法模型
x.fit <- decompose(x,type='addi')

# 可以查看季节指数、趋势效应、随机效应
x.fit$figure
x.fit$trend
x.fit$random

无季节效应差分运算

1
2
diff(x,d,k)
# 进行k次d步差分

差分

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
x <- ts()

# 画原时间序列
plot(x)
# 原时序图有线性趋势

# 1阶差分
x.dif<-diff(x)

# 1阶差分后的时序图\自相关图和偏相关图
plot(x.dif)
# 时序图在均值附近稳定波动

acf(x.dif)
pacf(x.dif)
# 观察截尾拖尾性质，判断拟合模型


ARIMA模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
arima(x,order=,include.mean=,method=,transform.par=,fixed=,seasonal=)
# order为非季节效应部分的模型阶数
# 是否需要拟合常数项，参数估计方法，是否需要人为干预参数估计，疏系数模型指定疏系数的位置
# seasonal =  list(order=c(P,D,Q),period=s)
# 加法：p=0,q=o;乘法:p/q不全为0
# 拟合
x.fit<-arima(x,orderc(0,1,1))
# 有季节因素的乘法模型
fit2 <-arima(x3,order =c(0,0,1),seasonal= list(order =c(0,1,1)),period =12),method='ML')#经过多次试探性降阶尝试，拟合AIC值较
# 残差白噪声检验
for (i in 1:2) print(Box.test(x.fit,lag=6*i))
# 白噪声p值大于0.05：显著成立才算拟合成功

# 预测
x,fore<-forecast(x.fit,h=10)

# 预测图

根据拟合结果写出拟合模型

1
2
x.fit
# 观察输出

疏系数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
arima(x,order=,include.mean=,method=,transform.pars=,fixed=)
-x:要进行模型拟合的序列名，
-order:
指定模型阶数。order=c(p,d,q),p为自回归阶数，d为差分阶数，q为移动平均阶数。
-include,mean:
include.mean=T,需要拟合常数项；include.mean=F,不需要拟合常数项，
-method:指定参数估计方法，
-transform,pars:指定参数估计是否由系统自动完成，
(1)transform.pars=T:系统默认设置是系统根据order选项设置的模型阶数自动完成参数估计。
(2)transform,pars=F:需要拟合疏系数模型，不能让系统根据模型的最高阶数自动完成所有参数的估计，我们需要进行人为干预，
-fixed：对疏系数模型指定疏系数的位置。

# 若偏自相关图，除了1阶和4阶偏自相关系数显著大于2倍标准差，其他阶数的偏自相关系数基本都在2倍标差范围内波动，疏系数AR(1,4)
# 疏系数模型
## ARIMA((1,4),1,0)
fit <- arima(x,order=c(4,1,0),transform.pars=F,fixed=c(NA,0,0,NA))
# transform.pars=F 人为干预
#fixed=c(NA,0,0,NA):只有1，4两个参数非零；2，3两个参数恒等于0
## # ARIMA(3,2,(1,2,6,7))
g.fit4<-arima(lng,order=c(3,2,7),transform.pars = F,fixed = c(NA,NA,NA,NA,NA,0,0,0,NA,NA))# 包含了p的三个系数

移动平均提取趋势

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 简单移动平均
# 简单移动平均就是将n个观测值的平均数作为第（n+1）/2个的拟合值。当n为偶数时，需进行二次移动平均。
m4 <- fliter(x/5,rep(1,5))


# 2*4复合移动平均

# 先4（周期长度）
m4 <- fliter(x/4,rep(1,4))
# 再2
m2_4 <- fliter(x4/2,rep(1,2),sides=1)
# sides = 1或者2，“1”表示单边卷积，“2”表示双边卷积
# 显示结果
dataframe(x,m2_4)

# 绘制移动效果
plot(x,lty=2)
lines(m2_4,col=2)

# 残差序列图
# 原序列-趋势效应 = 季节效应+随机波动
x_t <- x-m2_4
plot(x_t)

加法季节

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# x.t为剔除趋势效应的原序列
x.t<-

# 剔除缺失值，求总均值m
m<-mean(x_t,na.rm=T)

# 求每个季节的均值ms
ms <- 0
for (k in 1:4)ms[k]=mean(x_t[,k],na.rm=T)

# 季节指数
S<-ms-m

# 四个季节指数
plot(c(1:4),S,type='o')

# 随机效应
I<-x-m2_4-S
plot(I)

乘法模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# x.t为剔除趋势效应的原序列
x.t<-

# 剔除缺失值，求总均值m
m<-mean(x_t,na.rm=T)

# 求每个季节的均值ms
ms <- 0
for (k in 1:4)ms[k]=mean(x_t[,k],na.rm=T)

# 季节指数
S<-ms/m

# 12个季节指数
Month <- seq(1:12)
plot(Month,S,type='o')

# 随机效应
I<-x/m2_12/S
plot(I)

异方差模型

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
library(moments)
library(rugarch)
library(fGarch)
library(FinTS)
#LM检验(拉格朗日乘子检验)
for(i in 1:5) print(ArchTest(lnr_d,lag=i))
# LM检验显示ARCH模型的P值均小于0.05，显著成立，这说明残差平方序列具有短期相关性，可以用ARCH模型提取残差平方序列中蕴含的相关关系。
# 拟合arch(1)模型
spec2<-ugarchspec(mean.model = list(armaOrder=c(0,0),include.mean=F),variance.model = list(garchOrder=c(1,0),model="sGARCH"),distribution.model="norm")
fit4<-ugarchfit(spec2,data=lnr_d,method="ML")
fit4
fit4@fit$matcoef
plot(fit4,which=8)
plot(fit4,which=9)
plot(fit4,which=1,ylim=c(-300,300))
abline(h=c(-1.96*sd(lnr),1.96*sd(lnr)),col=1,lwd=2,lty=2)
fore <- ugarchforecast(fit4,5)

多元

例子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
library('TSA')
arimax(y,order=c(p,d,q),xreg=,xtransf=,transfer=)
#xreg:输入变量名（不需做函数转移）
#xtransf：输入变量名（需做函数转移）
#transfer：指定转移函数的模型阶数

y.fit=arimax(y,order=c(p,d,q),xreg=x)
#单变量不需要转移函数

y.fit=arimax(y,order=c(p,d,q),xreg=data.frame(x1,x2))
#两变量不需要转移函数

y.fit=arimax(y,order=c(p,d,q),xtransf=x,transfer=list(c(m,n)))
#单变量，ARMA(m,n)转移函数

y.fit=arimax(y,order=c(p,d,q),xtransf=data.frame(x1,x2),transfer=list(c(m1,n1),c(m2,n2)))
#双变量，分别ARMA(m,n)函数

y.fit=arimax(y,order=c(p,d,q),xreg=x1,xtransf=x2,transfer=list(c(m,n))
#双变量，一个需要转移，一个不需要转移。

x1<-zlag(x,k)
x1<-x1[-c(1:k)]
y1<-y[-c(1:k)]
#单变量，延迟k阶再ARMA(m,n)转移函数
-------------未完待续欢迎下次光临-------------
谢谢你想打赏的心~