#差分运算

#差分示例

#ARIMA模型

#Box-Jenkins建模流程

#案例

实际上我们经常会遇到一些非平稳时间序列，往往会呈现明显的趋势性或周期性，可以通过适当差分等手段，将它化为平稳时间序列，再采用ARMA模型建模

• 差分运算

差分运算公式

一阶差分 

二阶差分 

d 阶差分  ，其中 B 表示后移算子，

季节差分 

D阶季节差分  

差分示例

有时候差分还是不能够让序列平稳，那么我就可以使用对数转换等方式，然后再进行差分

• ARIMA模型

• 差分方式的选择

序列蕴含着显著的线性趋势，一阶差分就可以实现趋势平稳

序列蕴含着曲线趋势，通常低阶（二阶或三阶）差分可以提取出曲线趋势的影响，有时则需要对序列做非线性变换

对于蕴含着固定周期的序列进行步长为周期长度s的差分运算（季节差分），可以较好地提取周期信息

• Box-Jenkins建模流程

模型识别

自相关系数 偏

自相关系数

模型估计

最小化残差平方和

最小化AIC等等

模型验证

一些诊断方法会被用于检测模型的有效性

残差白噪声检验，模型系数显著性检验

用验证集数据来验证

模型预测

模型确定，用于预测

案例

> plot(tsdata)  #做时序图观察，是否为平稳的序列，其实可以直接adf.test() 检测，如果显著则满足平稳，如下图

> log_tsdata<-log(tsdata)  #因为数据变化幅度大，恐直接差分可能无法平稳，所以非线性变换，如进行log变换
> plot.ts(log_tsdata) #再绘制图形如下，发现变化幅度减小了

> log_tsdata_diff1=diff(log_tsdata,lag=12,differences=1)#lay=12表示差分周期（步长）12，differences=1表示一阶
> par(mfrow=c(3,1))
> plot.ts(log_tsdata_diff1)
> acf(log_tsdata_diff1)  #绘制 acf、pacf确定p、q的值
> pacf(log_tsdata_diff1)
> adfTest(log_tsdata_diff1,type='nc') #平稳性检验Title:Augmented Dickey-Fuller TestTest Results:PARAMETER:Lag Order: 1STATISTIC:Dickey-Fuller: -1.3524P VALUE:0.1828    #不显著，结果并不平稳，在进行处理一次Description:Wed Oct 18 09:47:53 2017 by user: Xu> log_tsdata_diff2=diff(log_tsdata_diff1,lag=1,differences=1)#步长为1的普通差分
> par(mfrow=c(3,1))
> plot.ts(log_tsdata_diff2)
> acf(log_tsdata_diff2)
> pacf(log_tsdata_diff2)
> adfTest(log_tsdata_diff2,type='nc')Title:Augmented Dickey-Fuller TestTest Results:PARAMETER:Lag Order: 1STATISTIC:Dickey-Fuller: -9.3424P VALUE:0.01   #显著表示平稳

 #建立模型并评价模型
> arima.mod=arima(log(tsdata), c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12)) #preiod周期为12#第一个c(0,1,1)表示一阶普通差分#第二个seasonal = list(order = c(0, 1, 1), period = 12)表示按照周期12的步长，进行季节性差分
> arima.modCall:
arima(x = log(tsdata), order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), period = 12))Coefficients:ma1     sma1-0.4018  -0.5569
s.e.   0.0896   0.0731sigma^2 estimated as 0.001348:  log likelihood = 244.7,  aic = -483.4> resid=arima.mod$residuals  #模型进行评价
> Box.test(resid)Box-Pierce testdata:  resid
X-squared = 0.03, df = 1, p-value = 0.8624> par(mfrow=c(1,1))
> qqnorm(resid)
> qqline(resid)

#预测数据
> library(forecast)
> arima.pred<-forecast(arima.mod,h=5*12) #预测5年的目标
> plot.forecast(arima.pred)
> ts.plot(tsdata,exp(arima.pred$mean), lty = c(1,3))
> AIC(arima(log(tsdata), c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12)))#通过AIC判断那个模型参数最好，可手动设置c(p,d,q)中的参数进行尝试，选取AIC最小的
[1] -483.3991
> AIC(arima(log(tsdata), c(0, 1, 2),seasonal = list(order = c(0, 1, 1), period = 12)))
[1] -481.6165

自动选择模型

> arima.new<-auto.arima(log(tsdata))#直接通过auto.arima来判断参数最合适
> arima.new
Series: log(tsdata)
ARIMA(0,1,1)(0,1,1)[12]                    Coefficients:ma1     sma1-0.4018  -0.5569
s.e.   0.0896   0.0731sigma^2 estimated as 0.001371:  log likelihood=244.7
AIC=-483.4   AICc=-483.21   BIC=-474.77