平稳时间序列分析
方法性工具
差分运算
- 一阶差分
- p阶差分
- k步差分
延迟算子
线性差分方程
ARMA模型
AR模型
-
AR模型平稳性判别
- 特征根判别
- 平稳域判别
-
平稳AR模型的统计性质
MA模型
-
MA模型的统计性质
- 均值
- 方差
- 自协方差q阶截尾
- 自相关系数q阶截尾
-
MA模型的可逆性
-
偏自行管系数截尾
ARMA模型
- 平稳条件 & 可逆条件
- 传递形式 & 逆转形式
- ARMA(p,q)模型的统计性质
平稳序列建模
建模步骤
- 求序列的样本自相关系数ACF和样本偏自相关系数PACF值
- 根据两个系数的性质,选择阶数适当的ARMA(p,q)模型进行拟合
- 估计模型中未知参数的值
- 检验模型有效性。如果拟合不符,返回步骤2重新选择模型拟合
- 模型优化。通过检验,返回步骤2,充分考虑各种可能,建立多个拟合模型,选择最优模型
- 利用拟合模型,预测序列将来走势
样本自相关系数 & 偏自相关系数
模型识别
参数估计
- 矩估计
- 极大似然估计
- 最小二乘估计
模型检验
- 模型显著性检验
- 参数显著性检验
模型优化
- 问题提出
- AIC准则
- SBC(BIC)准则
序列预测
线性预测函数
预测方差最小原则
线性最小方差与预测的性质
- 条件无偏最小方差估计值
- AR§序列预测
- MA(q)序列预测
- ARMA(p,q)序列预测
#AR模型平稳性判别
#1、arima.sim函数拟合
#arima.sim(n,list(ar=,ma=,order=),sd=)
#n:拟合序列长度
#list:指定具体的模型参数
#(1)拟合平稳AR(p)模型,要给出自回归系数,如果指定拟合的AR模型为非平稳模型,系统会报错
#(2)拟合MA(q)模型,要给出移动平均系数
#(3)拟合平稳ARMA(p,q)模型需要同时给出自回归系数和移动平均系数,如果指定拟合的ARMA模型为非平稳模型,系统会报错
#(4)拟合ARIMA(p,d,q)模型(第5章介绍),除了需要给出自回归系数和移动平均系数,还需要增加order选项。order=c(p,d,q)
# 其中,p为自回归阶数;d为差分阶数;q为移动平均阶数
#sd:指定序列的标准差,不特殊指定,系统默认 sd=1
#2、filter函数拟合
#filter(e,filter=,method=,circular=)
#e:随机波动序列的变量名
#filter:指定模型系数
#(1)AR(p)模型为filter=c(∅1,∅2,…,∅p)
#(2)MA(q)模型为filter=c(1,-θ1,-θ2,…,-θq)
#method:指定拟合的是AR模型还是MA模型
#(1)method=“recursive”为AR模型
#(2)method=“convolution”为MA模型
#circular:拟合MA模型时专用的一个选项,circular=T可以避免NA数据出现
#3.1
x1=arima.sim(n=100,list(ar=0.8))
x3=arima.sim(n=100,list(ar=c(1,-0.5)))
e=rnorm(100)
x2=filter(e,filter=-1.1,method = "recursive")
x4=filter(e,filter = c(1,0.5),method="recursive")
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
#layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
# layout() 输入一个矩阵,2行2列,
#第1个图占两个1,第2个图占2,第3个图占3
ts.plot(x1) #平稳
ts.plot(x2) #非平稳
ts.plot(x3) #平稳
ts.plot(x4) #非平稳
#3.5
x1=arima.sim(n=1000,list(ar=0.8))
x2=arima.sim(n=1000,list(ar=-0.8))
x3=arima.sim(n=1000,list(ar=c(1,-0.5)))
x4=arima.sim(n=1000,list(ar=c(-1,-0.5)))
#AR模型样本自相关图
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
acf(x1)
acf(x2)
acf(x3)
acf(x4)
#AR模型样本偏自相关图
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
pacf(x1)
pacf(x2)
pacf(x3)
pacf(x4)
#3.6
x1=arima.sim(n=1000,list(ma=-2))
x2=arima.sim(n=1000,list(ma=-0.5))
x3=arima.sim(n=1000,list(ma=c(-4/5,16/25)))
x4=arima.sim(n=1000,list(ma=c(-5/4,25/16)))
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
acf(x1)
acf(x2)
acf(x3)
acf(x4)
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
pacf(x1)
pacf(x2)
pacf(x3)
pacf(x4)
#3.8
x=arima.sim(n=1000,list(ar=0.5,ma=-0.8))
acf(x)
pacf(x)
#3.9
#读入数据,并绘制时序图
layout(matrix(c(1), 1, byrow = TRUE))
a=read.table("E:/data/file8.csv",sep=",",header=T)
x=ts(a$kilometer,start=1950)
plot(x)
#白噪声检验
for(i in 1:2) print(Box.test(x,type="Ljung-Box",lag=6*i))
#绘制自相关图和偏自相关图
acf(x)
pacf(x)
#3.10(有误)
#读入数据,并绘制时序图
overshort<-read.table("E:/data/file9.csv",sep=",",header=T)
overshort<-ts(overshort)
plot(overshort)
#白噪声检验
for(i in 1:2) print(Box.test(overshort,type="Ljung-Box",lag=6*i))
#绘制自相关图和偏自相关图
acf(overshort) #一阶截尾
pacf(overshort) #拖尾
#3.11
#读入数据,并绘制时序图
b=read.table("E:/data/file10.csv",sep=",",header=T)
dif_x=ts(diff(b$change_temp),start=1880)
plot(dif_x)
#白噪声检验
for(i in 1:2) print(Box.test(dif_x,type="Ljung-Box",lag=6*1))
#绘制自相关图和偏自相关图
acf(dif_x)
pacf(dif_x)
#auto.arima函数:先安装packages:zoo & forecast,然后用library调用程序包
#auto.arima(x,max.p=5,max.q=,ic=)
#x:需要定阶的序列名
#max.p:自相关系数最高阶数,不特殊指定的话,系统默认值为5
#max.q:移动平均系数最高阶数,不特殊指定的话,系统默认值为5
#ic:指定信息量准则,ic有"aicc","aic"和"bic"三个选项,系统默认AIC准则
library(zoo)
library(forecast)
#3.9 系统自动定阶
auto.arima(x)
#3.10 系统自动定阶
auto.arima(dif_x)
#3.11 系统自动定阶
auto.arima(dif_x)
#arima(x,order=,include.mean=,method=)
#x:要进行模型拟合的序列名
#order:指定模型阶数。order=c(p,d,q)
#(1)p为自回归阶数
#(2)d为差分阶数,本章不涉及缠粉问题,所以d=0
#(3)q为移动平均阶数
#include.mean:要不要包含常数项
#(1)include.mean=T,需要拟合常数项,这也是系统默认设置
#(2)如果不需要拟合常数项,需要特别指定include.mean=F
#method:指定参数估计方法
#(1)method="CSS-ML",默认的是条件最小二乘与极大似然估计混合方法
#(2)method="ML",极大似然估计
#(3)mrthod="CSS",条件最小二乘估计
#3.9续
a=read.table(file="E:/data/file8.csv",sep=",",header=T)
x=ts(a$kilometer,start=1950)
x.fit=arima(x,order=c(2,0,0),method="ML")
x.fit
#3.10续(有误)
overshort=read.table("E:/data/file9.csv",sep=",",header = T)
overshort=ts(overshort)
overshort.fit=arima(overshort,order=c(0,0,1))
overshort.fit
#3.11续
b=read.table("E:/data/file10.csv",sep=",",header = T)
dif_x=ts(diff(b$change_temp),start = 1880)
dif_x.fit=arima(dif_x,order=c(1,0,1))
dif_x.fit
#3.9续
a=read.table(file="E:/data/file8.csv",sep=",",header = T)
x=ts(a$kilometer,start = 1950)
x.fit=arima(x,order=c(2,0,0),method = "ML")
for(i in 1:2) print(Box.test(x.fit$residuals,lag=6*i))
#3.10续(有误)
overshort=read.table("E:/data/file9.csv",sep=",",header = T)
overshort=ts(overshort)
overshort.fit=arima(overshort,order=c(0,0,1))
for(i in 1:2) print(Box.test(overshort.fit$residual,lag=6*i))
#3.11续
b=read.table("E:/data/file10.csv",sep=",",header = T)
dif_x=ts(diff(b$change_temp),start = 1880)
dif_x.fit=arima(dif_x,order=c(1,0,1),method="CSS")
for(i in 1:2) print(Box.test(dif_x.fit$residual,lag=6*i))
#pt(t,df=,lower.tail=)
#t:t统计量的值
#df:自由度
#lower.tail:确定计算概率的方向
#(1)lower.tail=T,计算Pr(X≤x)。对于参数显著性检验,如果参数估计值为负,选择lower.tail=T
#(2)lower.tail=F,计算Pr(X>x)。对于参数显著性检验,如果参数估计值为正,选择lower.tail=F
#3.9续
a=read.table(file="E:/data/file8.csv",sep=",",header = T)
x=ts(a$kilometer,start = 1950)
x.fit=arima(x,order=c(2,0,0),method="ML")
x.fit
#3.15(有误)
#读入数据,绘制时序图
x=read.table(file="E:/data/file11.csv",sep=",",header = T)
x=ts(x)
plot(x)
#序列白噪声检验
for(i in 1:2) print(Box.test(x,lag=6*i))
#绘制自相关图和偏自相关图
acf(x)
pacf(x)
#拟合MA(2)模型
x.fit1=arima(x,order=c(0,0,2))
x.fit1
#MA(2)模型显著性检验
for(i in 1:2) print(Box.test(x,fit1$residual,lag=6*1))
#拟合AR(1)模型
x.fit2=arima(x,order=c(1,0,0))
x.fit2
#AR(1)模型显著性检验
for(i in 1:2) print(Box.test(x.fit2$residual,lag=6*i))
#forecast(objuct,h=,level=)
#object:拟合信息文件名
#h:预测期数
#level:置信区间的置信水平。不特殊指定的话,系统会自动给出置信水平分别为80%和95%的双层置信区间
#3.9续
a=read.table(file="E:/data/file8.csv",sep=",",header = T)
x=ts(a$kilometer,start = 1950)
x.fit=arima(x,order=c(2,0,0))
x.fore=forecast(x.fit,h=5)
x.fore
#系统默认输出预测图
plot(x.fore)
#个性化输出预测图
L1=x.fore$fitted-1.96*sqrt(x.fit$sigma2)
U1=x.fore$fitted+1.96*sqrt(x.fit$sigma2)
L2=ts(x.fore$lower[,2],start = 2009)
U2=ts(x.fore$upper[,2],start = 2009)
c1=min(x,L1,L2)
c2=max(x,L2,U2)
plot(x,type="p",pch=8,xlim=c(1950,2013),ylim=c(c1,c2))
lines(x.fore$fitted,col=2,lwd=2)
lines(x.fore$mean,col=2,lwd=2)
lines(L1,col=4,lty=2)
lines(U1,col=4,lty=2)
lines(L1,col=4,lty=2)
lines(L2,col=4,lty=2)
lines(U2,col=4,lty=2)