时间序列分析

借助实例，理解R代码，形成自己的代码库，方便自己使用。

---

前言

一般地，需要将数据转化为时间序列数据。

library(tseries) #导入包

#构造时间序列
yield<- ts(a$yield,start=1964)
people<- ts(a,frequency = 4,start=c(1971,9))

---

一、平稳时间序列

平稳性检验

library(tseries) #导入包
x<- ts(x, start=1964)

# 时序图检验
plot(x)
# 自相关图检验
acf(x, lag.max=n) #lag.max最大的序列相关

# 做纯随机检验
# Box.test(x, type= , lag=)
# type="Box-Pierce #Q统计量, 默认
# type="Ljung-Box" #LB统计量
white_noise<-rnorm(100)
white_noise <-ts(white_noise)
Box.test(white_noise,lag=6)
for(i in 1:3) print(Box.test(white_noise,lag=3*i))

arima.sim函数

x1<-arima.sim(n=100,list(ar=0.8))
ts.plot(x1)#arima函数
x3<-arima.sim(n=100,list(ar=c(1,-0.5)))
ts.plot(x3)

e<-rnorm(100)
x2 <-filter(e, filter=-1.1, method='recursive')
ts.plot(x2)#fileter函数
x4 <-filter(e, filter=c(1,0.5), method='recursive')
ts.plot(x4)

#随机游走
x<-arima.sim(n=1000,list(order=c(0,1,0)),sd=10)
ts.plot(x)

差分方程求解

a <- c(a0,a1,z2, . ,ap)#输入差分方程a0*zt+a1*zt-1+...+ap*zt-p=0的系数
Φ <- c(-1,Φ1, . ,Φp)#输入系数时需要取相反数

polyroot(a)#求解方程z0+z1*x+z2*x^2+...+zp*x^p=0
eigen <- 1/ployroot(a)#取倒数得到特征根

ARMA模型

参数估计

# 自动定阶
library(zoo)
library(forecast)

auto.arima(x) #x是一个时间序列数据

# 手动设置
acf(x) #MA有限，则拖尾
pacf(x) #p阶截尾，则AR为p阶数
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))

预测

x.fore<-forecast(x.fit,h=5)
x.fore
plot(x.fore)

#个性化输出预测图
L1<-x.fore$fitted-1.96*sqrt(x.fit$sigma2)#95%置信区间下限,其中x.fore$fitted往期拟合的数据值
U1<-x.fore$fitted+1.96*sqrt(x.fit$sigma2)
L2<-ts(x.fore$lower[,2],start=2009)#95%置信区间下限
U2<-ts(x.fore$upper[,2],start=2009)
c1=min(x,L1,L2);c2=max(x,L1,L2)#确定众坐标取值范围
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(L2,col=4,lty=2)
lines(U2,col=4,lty=2)

二、非平稳序列确定性分析

趋势拟合法

lm(x~t) #线性拟合
nls(x~f(t), data=, start=) #曲线拟合

平滑法

#移动平均法
library(TTR)
x<-ts(a$temp,start = 1949)
x.ma<-SMA(x,n=5)
plot(x,type='o')
lines(x.ma,col=2,lwd=2)

#指数平滑法
HoltWinters(x, alpha=NULL, beta=NULL, gamma=NULL, seasonal="additive" )

x<-ts(a$output,start = c(1962,1),frequency = 12)
x.fit<-HoltWinters(x)
x.fit
# 预测
library(forecast)
x.fore<-forecast(x.fit,h=10)
x.fore
plot(x.fore)

季节效应分析

x<-ts(a$output,start = c(1962,1),frequency = 12)
x.fit<-decompose(x,type='mult')
x.fit$figure
plot(x.fit$figure,type='o')
x.fit$trend
plot(x.fit$trend)
x.fit$random
plot(x.fit$random)
# 综合图
plot(x.fit)

三、非平稳序列随机性分析

差分

为了使序列平稳，可以进行差分。原则只要进行无限的差分，序列总能平稳，但过差分会使得残差序列的方差变大。

# 观察时序图，确定差分的阶数和周期
# 一阶差分
diff(x)
# 12步差分
diff(x, 12)

# 使用ARIMA模型进行拟合
auto.arima(x) #自动定阶
arima(x, order=c(0,1,1))
#预测与ARMA一致

疏系数模型

序列可能不存在一阶自相关，但是存在二阶自相关，这时需要疏系数模型。

arima(x, order=c(4,1,0), transform.pars=F, fixed=c(NA,0,0,NA))

季节模型

#加法季节模型
arima(x, order=c(4,1,0), seasonal=list(order=c(0,1,0), period=4), transform.par=F, fixed=c(NA,0,0,NA)) 
#乘法季节模型
arima(x, order=c(4,1,0), seasonal=list(order=c(0,1,0), period=12)) 

残差自回归模型

#拟合趋势
x.fit <- lm(x~f(t))
x.fit$fitted.value #趋势预测图
#得到残差
res <- x.fit$residual
auto.arima(res)

异方差

lnx<-log(x) #变换函数
plot(lnx)
dif.lnx<-diff(lnx)
plot(dif.lnx)
for(i in 1:2)print(Box.test(dif.lnx,lag=6*i))

ARCH(集群效应)

观察序列

x<-ts(a$exchange_rates, start=c(1979,12),frequency = 12)
plot(x^2)

拟合ARCH（3）

x.fit<-garch(x, order = c(0,3))
summary(x.fit)
x.pred<-predict(x.fit)
plot(x.pred)#绘制条件异方差模型拟合的95%置信区间
#条件异方差置信区间和方差齐性的置信区间
plot(x)
lines(x.pred[,1],col=2)
lines(x.pred[,2],col=2)
abline(h=1.96*sd(x),col=4,lty=2)
abline(h=-1.96*sd(x),col=4,lty=2)

AR-GARCH模型

x<-ts(a$exchange_rates,start=c(1997,12,31),frequency = 365)
x.fit<-arima(x, order = c(0,1,1))
# 拟合残差
for(i in 1:6) print(Box.test(x.fit$residual^2,type="Ljung-Box",lag=i))
#拟合GARCH(1,1)模型
r.fit<-garch(x.fit$residuals,order=c(1,1))
summary(r.fit)
#绘制波动置信区间
r.pred<-predict(r.fit)
plot(r.pred)

四、多元时间序列

多元时间序列回归

y.fit <- arimax(y, order=c(p,0,q), xreg=x) # 只有y需要进行变换
y.fit <- arimax(y, order=c(p,0,q), xreg=data.frame(x1, x2))
# y, x1, x2 都需要进行变换
y.fit <- arimax(y, order=c(p,0,q), xtransf=data.frame(x1,x2), transfer=list(c(m1,n1), c(m2,n2))

平稳性检验

library(fBasics)
library(fUnitRoots)
#DF检验
adfTest(x,lag=1,type = "nc") #无常值平稳
adfTest(x,lag=1,type = "c")	#有常值平稳
adfTest(x,lag=1,type = "ct") #有趋势平稳
#ADF检验1, 2, 3阶是否平稳
for(i in 1:3) print(adfTest(x,lag=i,type="nc"))
for(i in 1:3) print(adfTest(x,lag=i,type="c"))
for(i in 1:3) print(adfTest(x,lag=i,type="ct"))

# 一阶差分平稳
dx <- diff(x)

协整

检验回归残差序列是否平稳

y.fit <- lm(y~x)
r <- ts(y.fit$residual)
for(i in 1:3) print(adfTest(r, lag=i, type="nc"))

误差修正模型

y.fit <- lm(y~x)
ECM <- y.fit$residual[1:length(x)-1]
dify <- lm(diff(y)~0+diff(x)+ECM)
summary(dify)