非平穩時間序列及建模

#差分運算

#差分示例

#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 Test

Test Results:
  PARAMETER:
    Lag Order: 1
  STATISTIC:
    Dickey-Fuller: -1.3524
  P 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 Test

Test Results:
  PARAMETER:
    Lag Order: 1
  STATISTIC:
    Dickey-Fuller: -9.3424
  P 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.mod

Call:
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.0731

sigma^2 estimated as 0.001348:  log likelihood = 244.7,  aic = -483.4


> resid=arima.mod$residuals  #模型進行評價
> Box.test(resid)

	Box-Pierce test

data:  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.0731

sigma^2 estimated as 0.001371:  log likelihood=244.7
AIC=-483.4   AICc=-483.21   BIC=-474.77