R时间序列分析

R时间序列分析   为什么定阶数,如何定,如何判断
R时间序列分析工具
xts包  xts(x=NUll,order.by=index(x),…)  coredata()   xts数据子集
OHLC数据格式
quantmod包
TTR包
自回归模型(AR)  跟以前时刻有关和当前随机游动有关
AR(p)的性质
平稳性要求：
AR(2)特征根的模都小于1
AR(p)的定阶 赤池信息准则（AIC） 施瓦茨—贝叶斯信息准则（BIC） 
AIC BIC最小值的地方定阶
利用偏自相关函数（PACF）与上一刻时间的关系  与前几个时刻相关的显著性 越大越显著
白噪声时是严平稳序列
平稳时间序列
# AR models 自回归模型AR(p) p阶
da=read.table("q-gnp4710.txt",header=T)
head(da)
G=da$VALUE
LG=log(G)  #取对数
gnp=diff(LG)  #LG的差分,返回适当的滞后，迭代的差异
dim(da)  #重新设置或设置对象的维数
tdx=c(1:253)/4+1947 # create the time index
par(mfcol=c(2,1))
plot(tdx,G,xlab='year',ylab='GNP',type='l')
plot(tdx[2:253],gnp,type='l',xlab='year',ylab='growth') 
acf(gnp,lag=12)  #计算互相关或两个单变量序列的协方差
pacf(gnp,lag=12) # compute PACF
m1=arima(gnp,order=c(3,0,0)) #拟合ARIMA模型的单变量时间序列  order=c(3,0,0)第一个是AR阶数第二个是差分阶数第三个是MA模型阶数
m1
tsdiag(m1,gof=12)  # model checking discussed later绘制时间序列
p1=c(1,-m1$coef[1:3]) # set-up the polynomial设置多项式
r1=polyroot(p1) # solve the polynomial equation求解多项式方程
r1
Mod(r1)  #绝对值 (模)
k=2*pi/acos(1.616116/1.832674) # compute length of the period 计算周期长度
k
Box.test(m1$residuals,lag=4,type='Ljung')  #模型检验
1-pchisq(0.74926,1)
forecast(m1,1) #预测第一步
forecast(m1,2) #预测第二步
滑动平均模型(MA) 与上几个随机游动有关
da=read.table("m-ibm3dx2608.txt",header=T)
head(da)
ew=da$ewrtn
plot(da$data,ew,type="1")
#定阶  取不为0的阶不要太高(复杂)不要太低(不准)非0说明相关
acf(ew,lag=20)
#模型拟合
m1=arima(ew,order=c(0,0,9)) 
m1
m1=arima(ew,order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))#只拟合相关系数不为0的项  aic越小越好
m1
sqrt(0.005097)  #计算随机游动标准差
#模型检验
tsdiag(m1,gof=12)  #检测残差是否为白噪声，残差为白噪声好
Box.test(m1$residuals,lag=12,type='Ljung')  # model checking
pv=1-pchisq(17.6,9)  # compute p-value after adjusting the d.f.
pv
#预测
m1=arima(ew[1:986],order=c(0,0,9),fixed=c(NA,0,NA,0,0,0,0,0,NA,NA))#只拟合相关系数不为0的项  aic越小越好
m1
predict(m1,10)
pre=c(predict(m1,10)$pred)
data.frame(ew[987:996],pre,ew[987:996]-pre)  #预测值与真实值比较
自回归滑动平均模型(ARMA)
AR(p)模型与MA(q)模型的结合——ARMA(p,q)
ARMA(p,q)定阶
推广的自相关函数（EACF） 非零
#构建ARMA模型
da=read.table("m-ibm3dx2608.txt",header=T)
head(da)
ew=da$ewrtn
acf(ew)
pacf(ew)
#根据信息准则定阶
library("forecast")
auto.arima(ew)
#根据eacf定阶
eacf(ew)
#模型拟合
m2=arima(ew,order=c(2,0,2))
m2
m3=arima(ew,order=c(0,0,1))
m3
#模型检验
tsdiag(m2,gof=12)  #检测残差是否为白噪声
Box.test(m2$residuals,lag=12,type='Ljung')  # model checking
tsdiag(m3,gof=12)  #检测残差是否为白噪声
Box.test(m3$residuals,lag=12,type='Ljung')  # model checking
#预测
predict(m2,1)
predict(m3,1)
平稳性检验
ADF单位根检验library("fUnitRoots")
adfTest()
时间序列  确定性趋势成分+平稳的零均值误差成分
确定性因素分解
   长期趋势变动
 循环变动
 季节性变动
 随机波动
确定性时间序列分析
– 平滑预测法（平稳时间序列）
   简单移动平均法SMA
   k期移动平均法
   k期中心移动平均
– 指数平滑法
– 趋势预测法
   线性拟合
   曲线拟合
– 分解分析法
加权移动平滑（WMA）
– k期移动平均
– 指数加权移动平滑（EWMA）
简单移动平均法SMA
sma.cal <- function( ts ) {
n <- length(ts)
sma <- c()
for( t in 1 : n ) sma[t]= mean(ts[1:t])
return(sma)
}
加权移动平均法(时间越近,权重越高)
wma.cal <- function( ts, weight ){
n <- length(ts) #ts时间序列
k <- length(weight)  #weight权重数据
wma <- c()
for( t in (k+1):n ) wma[t]= sum( ts[(t-k+1):t] * weight )
return(wma)
}
k期移动平均法(简单移动平均法的优化,取最近K期的均值)
kwma.cal <- function( ts, k=20 ) {
n <- length(ts)
kwma <- c()
for( t in (k+1): n ) kwma[t]= mean(ts[(t-k+1):t])
return(kwma)
}
指数平滑法  平滑系数α的选择(时间越远权重小)
ewma.cal <- function( ts, a=0.8 ){
n <- length(ts)
ewma <- c(ts[1])
for( t in 2:n ) ewma[t]= a*ts[t] + (1-a)*ewma[t-1]
return(ewma)
}
aa=c(1:10)
n=length(aa)
a=0.8
pred=0
pred[1:3]<-aa[1:3]
for(i in 4:length(aa)){
pred[i]<-a*aa[i-1]+(1-a)*pred[i-1] #a是平滑指数，可人为控制
}#pred为预测值，其中前三个数据为真实数据
p=a*aa[n]+(1-a)*pred[n]

季节均值法
余弦趋势法
季节指数水平法
季节指数趋势法
全年比率平均法
步骤：
 1. 求各年平均数
 2. 求各月对该年平均数比率
 3. 计算季节比率
 4. 预测

非平稳时间序列
平稳性检验  ADF单位根检验 adfTest()
library(fUnitRoots)
adfTest(oil.price,lags=12,type="c")  #得出p值 取到这个样本的概率是多少 P值即概率，反映某一事件发生的可能性大小
参数lags(滞后项)的设置——根据AR模型的阶数估算而确定
data("oil.price") #TSA包里有
ar(oil.price,method = 'mle')$order  #得出阶数
识别趋势非平稳和差分非平稳
趋势非平稳
差分非平稳
library(fUnitRoots)
da=read.table("d-sp55008.txt",header=T)
head(da)
sp5=log(da[,7])
m2=ar(diff(sp5),method='mle') # Based on AIC
m2$order
adfTest(sp5,lags=2,type=("ct"))
adfTest(sp5,lags=15,type=("ct")) # Based on PACF
非平稳AR(1)模型  系数<1
data(explode.s)
plot(explode.s,ylab=expression(y[t]),type="o")
acf(explode.s)
pacf(explode.s)
对于随机游动 不易过度差分
一次差分后 是一个白噪声过程  在这里就可以停止了不易过度差分
两次差分后  MA(1)模型   过差分效果不好
ARIMA(p,d,q)模型
如果一个时间序列{