#加載數據

x=read.table(file.choose())

#生成時間序列對象

xtimeseries

#畫時間序列圖

plot.ts(xtimeseries)

#增加線性擬合曲線

abline(lm(xtimeseries~time(xtimeseries)))

1、分解時間序列

分解一個時間序列意味著把它拆分成構成元件，一般序列包含一個趨勢部分、一個不規則部分，如果是一個季節性時間序列，則還有一個季節性部分。

分解非季節性數據(趨勢部分+不規則部分)

估計趨勢部分:最常用的方法便是平滑法，比如計算時間序列的簡單移動平均。為了更加準確地估計這個趨勢部分，我們也許應該嘗試下更大的跨度進行平滑。

在R的“TTR”包中的SMA()函數可以用簡單的移動平均來平滑時間序列數據

# 使用SMA()函數時，你需要通過參數“n”指定來簡單移動平均的跨度

#這僅僅是一個代碼例子，該時間序列是有季節變化的

library(TTR)

xtimeseriesSMA

plot.ts(xtimeseriesSMA)

分解季節性數據(季節部分+趨勢部分+不規則部分)

#對于可以使用相加模型進行描述的時間序列中的趨勢部分和季節性部分，我們可以使用R中“decompose()”函數來估計。這個函數可以估計出時間序列中趨勢的、季節性的和不規則的部分，而此時間序列須是可以用相加模型描述的。

xtimeseriescomponents

plot(xtimeseriescomponents)

#從原始時間序列中去除去季節部分。

xtimeseriesseasonallyadjusted

# 我們可以使用“plot()”畫出季節性修正時間序列，代碼如下：

plot(xtimeseriesseasonallyadjusted)

如果你有一個可用相加模型描述的，并且處于恒定水平和沒有季節性變動的時間序列，你可以使用簡單指數平滑法對其進行短期預測。

簡單指數平滑法提供了一種方法估計當前時間點上的水平。為了準確的估計當前時間的水平，我們使用alpha參數來控制平滑。Alpha的取值在0到1之間。當alpha越接近0的時候，臨近預測的觀測值在預測中的權重就越小。

2、使用指數平滑法進行預測

指數平滑法可以用于時間序列數據的短期預測。

簡單指數平滑法(非季節性變動和處于恒定水平、沒明顯趨勢)

如果你有一個可用相加模型描述的，并且處于恒定水平和沒有季節性變動的時間序列，你可以使用簡單指數平滑法對其進行短期預測。

# 比如，使用簡單指數平滑發對倫敦每年下雨量進行預測，HoltWinters()函數中設定參數beta=FALSE和gamma=FALSE，代碼如下：

rainseriesforecasts

#預測的結果值

rainseriesforecasts$fitted

#樣本內預測誤差的誤差平方之和

rainseriesforecasts$SSE

使用R中的“forecast”包中的“forecast.HoltWinters()”函數進行更遠時間點上的預測。

library(forecast)

#h參數為預測之后的多少個時間窗口

rainseriescasts2

#繪畫預測結果

plot.forecast(rainseriescasts2)

霍爾特指數平滑法(有明顯趨勢，非季節的相加模型)

如果你的時間序列可以被描述為一個增長或降低趨勢的、沒有季節性的相加模型，你可以使用霍爾特指數平滑法對其進行短期預測。

Holt指數平滑法估計當前時間點的水平和斜率。其平滑化是由兩個參數控制的，alpha，用于估計當前時間點的水平，beta，用于估計當前時間點趨勢部分的斜率。

正如簡單指數平滑法一樣，alpha和beta參數都介于0到1之間，并且當參數越接近0，大多數近期的觀測則將占據預測更小的權重。

# 一個可能可以用相加模型描述的有趨勢的、無季節性的時間序列案例就是這1866年到1911年每年女人們裙子的直徑。這個數據可以從該文件獲得(http://robjhyndman.com/tsdldata/roberts/skirts.dat)(初始數據來自Hipel and McLeod, 1994)

skirts

skirtsseries

plot.ts(skirtsseries)

# 例如，為了使用Holt指數平滑法修正一個裙邊直徑的預測模型，我們輸入代碼：

skirtsseriesforecasts

skirtsseriesforecasts

# 這里的alpha預測值為0.84，beta預測值為1.00。這都是非常高的值，告訴我們無論是水平上，還是趨勢的斜率，當前值大部分都基于時間序列上最近的觀測值。

#這樣的直觀感覺很好，因為其時間序列上的水平和斜率在整個時間段發生了巨大的變化。預測樣本內誤差的誤差平方和是16954。

# 我們可以用黑色線條畫出原始時間序列分布，用紅色線條畫出頂部的預測值，代碼如下：

plot(skirtsseriesforecasts)

# 從該圖我們可以看到樣本內預測非常接近觀測值，盡管他們對觀測值來說有一點點延遲。

# 如果你想的話，你可以通過HoltWinters()函數中的“l.start”和“b.start”參數去指定水平和趨勢的斜率的初始值。常見的設定水平初始值是讓其等于時間序列的第一個值(在裙子數據中是608)，而斜率的初始值則是其第二值減去第一個值(在裙子數據中是9)。

# 例如，為了使用Holt指數平滑法找到一個在裙邊直徑數據中合適的預測模型，我們設定其水平初始值為608，趨勢部分的斜率初始值為9，代碼如下：

HoltWinters(skirtsseries,gamma=FALSE,l.start = skirtsseries[1],b.start = skirtsseries[2]-skirtsseries[1])

## Holt-Winters exponential smoothing with trend and without seasonal component.

##

## Call:

## HoltWinters(x = skirtsseries, gamma = FALSE, l.start = skirtsseries[1], b.start = skirtsseries[2] - skirtsseries[1])

##

## Smoothing parameters:

## alpha: 0.8346775

## beta : 1

## gamma: FALSE

##

## Coefficients:

## [,1]

## a 529.278637

## b 5.670129

# 正如我們的簡單指數平滑法一樣，我們可以使用“forecast”包中的forecast.HoltWinters()函數預測未來時間而無需覆蓋原始序列。

# 例如，我們的現在有的1866年到1911年的裙邊直徑時間序列數據，因此我們可以預測1912年到1930年(19個點或者更多)，并且畫出他們，代碼如下：

skirtsseriescasts2

plot.forecast(skirtsseriescasts2)

Holt-Winters指數平滑法(增長或降低趨勢、季節性變化、相加模型的時間序列)

如果你有一個增長或降低趨勢并存在季節性可被描述成為相加模型的時間序列，你可以使用霍爾特-溫特指數平滑法對其進行短期預測。

Holt-Winters指數平滑法估計當前時間點的水平，斜率和季節性部分。平滑化依靠三個參數來控制：alpha，beta和gamma，分別對應當前時間點上的水平，趨勢部分的斜率和季節性部分。

參數alpha，beta和gamma的取值都在0和1之間，并且當其取值越接近0意味著對未來的預測值而言最近的觀測值占據相對較小的權重。

#取對數可以減少極值帶來的影響，消除方差不齊。

logxtimeseries

xtimeseriesforecasts

xtimeseriesforecasts

#用黑線畫出原始數據的時間曲線圖，用紅線在上面畫出預測值的時間曲線圖：

plot(logxtimeseries)

plot(xtimeseriesforecasts)

#如果alpha相對較低，說明當前時間點估計得水平是基于最近觀測和歷史觀測值。beta的估計值是0.00，表明估計出來的趨勢部分的斜率在整個時間序列上是不變的，并且應該是等于其初始值。這是很直觀的感覺，水平改變非常多，但是趨勢部分的斜率b卻仍然是大致相同的。與此相反的，gamma的值(0.96)則很高，表明當前時間點的季節性部分的估計僅僅基于最近的觀測值。

# 我們可以從途中看出Holt-Winters指數平滑法是非常成功得預測了季節峰值，其峰值大約發生在每年的12月份。

# 為了預測非原始時間序列的未來一段時間，我們使用“forecast”包中的“forecast.HoltWinters()”函數。

# 例如，禮物消費原始數據是2013年7月到2016年8月。如果我們想預測2016年9月到2017年6月(12個月或者更多)，并且畫出預測，代碼如下：

xtimeseriesforecasts2

plot.forecast(xtimeseriesforecasts2)

# 我們可以通過畫相關圖和進行Ljung-Box檢驗來檢查樣本內預測誤差在延遲1-20階時否是非零自相關的，并以此確定預測模型是否可以再被優化

acf(xtimeseriesforecasts2$residuals,lag.max=20)

Box.test(xtimeseriesforecasts2$residuals,lag=20,type="Ljung-Box")

# Box-Ljung test

#

$data: xtimeseriesforecasts2$residuals

#X-squared = 16.4789, df = 20, p-value = 0.6865

# 這個樣本內預測誤差的相關圖并沒有在延遲1-20階內自相關系數超過置信界限的。而且，Ljung-Box檢驗的P值是0.68(>0.05)，意味著是不足以證明延遲1-20階是非零自相關的。

ps:自相關值超出95%的顯著邊界有可能是偶然的。

# 我們可以在整個時間段內檢驗預測誤差是否是方差不變，并且服從零均值正態分布。方法是畫出預測誤差的時間曲線圖和直方圖(并覆蓋上正太曲線)：

plot.ts(xtimeseriesforecasts2$residuals) # make a time plot

plotForecastErrors(xtimeseriesforecasts2$residuals) # make a histogram

ARIMA模型(AR代表自回歸，MA代表移動平均)

指數平滑法對于預測來說是非常有幫助的，而且它對時間序列上面連續的值之間相關性沒有要求。但是，如果你想使用指數平滑法計算出預測區間，那么預測誤差必須是不相關的，而且必須是服從零均值、方差不變的正態分布。

即使指數平滑法對時間序列連續數值之間相關性沒有要求，在某種情況下，我們可以通過考慮數據之間的相關性來創建更好的預測模型。自回歸移動平均模型(ARIMA)包含一個確定(explicit)的統計模型用于處理時間序列的不規則部分，它也允許不規則部分可以自相關。

時間序列的差分

ARIMA模型為平穩時間序列定義的。因此，如果你從一個非平穩的時間序列開始，首先你就需要做時間序列差分直到你得到一個平穩時間序列。如果你必須對時間序列做d階差分才能得到一個平穩序列，那么你就使用ARIMA(p,d,q)模型，其中d是差分的階數。

在R中你可以使用diff()函數作時間序列的差分。

# 我們可以通過鍵入下面的代碼來得到時間序列(數據存于“xtimeseries”)的一階差分，并畫出差分序列的圖：

xseriesdiff1

plot.ts(xseriesdiff1)

# 一階差分時間序列結果(上圖)的均值看起來平穩。因此，不需要再次做差分。如果不平穩繼續做差分。

#xseriesdiff2

#plot.ts(xseriesdiff2)

例如：

install.packages("tseries")

library(tseries)

data(AirPassengers)

adf.test(diff(log(AirPassengers)), alternative="stationary", k=0)

Augmented Dickey-Fuller Test

data: diff(log(AirPassengers))

Dickey-Fuller = -9.6003, Lag order = 0, p-value = 0.01

alternative hypothesis: stationary

#如果是這樣表示，這個序列是足夠平穩做任何時間序列模型。

# 對于平穩性正式的檢驗稱作“單位根測試”，可以在fUnitRoots包中得到

如果你的時間序列是平穩的，或者你通過做n次差分轉化為一個平穩時間序列，接下來就是要選擇合適的ARIMA模型，這意味著需要尋找ARIMA(p,d,q)中合適的p值和q值。(通過上面差分得出d=1)為了得到這些，通常需要檢查平穩時間序列的(自)相關圖和偏相關圖。

我們使用R中的“acf()”和“pacf”函數來分別(自)相關圖和偏相關圖。在“acf()”和“pacf設定“plot=FALSE”來得到自相關和偏相關的真實值。

acf(xseriesdiff1,lag.max=20) #plot a correlogram

acf(xseriesdiff1,lag.max=20,plot=FALSE) # get the autocorrelation values

# 我們從上面相關圖中可以看出在滯后1、6、12階(lag1)的自相關值(-0.360)超出了置信邊界。

# 要畫出禮物消費在滯后1-20階(lags 1-20)一階差分時間序列的偏相關圖，并得到偏相關的值，我們使用“pcf()”函數，輸入：

pacf(xseriesdiff1,lag.max=20) #plot a partial correlogram

# 偏相關圖顯示在滯后1，2和5階(lags 1,2,3)時的偏自相關系數超出了置信邊界，為負值。

確定pq

首先判斷acf圖和pacf圖是否平穩，加入假如非平穩那么需要差分，如果一階差分后仍非平穩，則需要二階差分，等等。在確定差分平穩后，需要判斷p和q，

這里定階方法有很多，因為p和q的確定也很復雜，不是一下子就可以確定的。主要有這么幾種(1)觀察法，直接觀察，如果acf在q+1階突然截斷，在q處截尾，則為ma(q)序列，同理，pacf 在p處截尾則為ar(p)序列，否則為arma(p,q)序列，二者結合進一步判斷

(2)參數檢驗，利用數理統計檢驗高階模型的新增加的參數是否近似為零，檢驗模型殘差的相關 特性等。

(3)信息準則，確定一個與模型階數有關的準則，如AIC、BIC等，既考慮擬合效果接近程度， 又考慮參數個數。實際中往往多種方法綜合應用，選擇最合適的p,d,q.

ARMA(p，q)模型的ACF與PACF理論模式

模型 ACF PACF

AR(p) 衰減趨于零(幾何型或振蕩型) p階后截尾(可以考慮鄰項觀測值上具有長期相關性)

MA(q) q階后截尾 衰減趨于零(幾何型或振蕩型)(移動平均模型鄰項觀察值之間短期相關的特征)

ARMA(p，q) q階后衰減趨于零(幾何型或振蕩型) p階后衰減趨于零(幾何型或振蕩型)

auto.arima() 函數可以用來發現合適的ARIMA模型，例如輸入 “library(forecast)”, 然后 “auto.arima(kings)”. 那么輸出說明合適的模型是ARIMA(0,1,1).

既然對于英國國王去世年齡的時間序列，ARIMA(0,1)被認為是最合適的模型，那么原始的時間序列可以使用ARIMA(0,1,1) (p=0,d=1,q=1，這里d是差分階層所需要的)來建模。

#auto.arima(logxtimeseries,ic='bic')如果你使用 “bic” 標準, 這里對參數個數要求非常嚴格, 使用簡單的原則(p+q最小)，模型在這里是同樣優先的選擇模型。

auto.arima(logxtimeseries,trace=T) #Best model: ARIMA(0,1,1)(0,1,1)[12]

xtimeseriesarima

你可以使用forecast.Arima() 中“level”參數來確定預測區間的置信水平。例如，為了得到99.5%的預測區間，我們輸入： “forecast.Arima(xtimeseriesarima, h=5, level=c(99.5))”.

然后我們可以使用ARIMA模型來預測時間序列未來的值，使用R中forecast包的“forecast.Arima()” 函數。例如，為了預測接下來5個英國國王的去世年齡，我們輸入：

library(forecast)

xtimeseriesforecasts

xtimeseriesforecasts

plot.forecast(xtimeseriesforecasts)

模型產生的預測誤差檢驗：

過畫相關圖和進行Ljung-Box檢驗來檢查樣本內預測誤差在延遲1-20階時否是非零自相關

acf(xtimeseriesforecasts$residuals,lag.max=20)

Box.test(xtimeseriesforecasts$residuals,lag=20,type="Ljung-Box") # p=0.4835>0.05

#Ps: p值大，說明為純隨機序列。p值小，非純隨機序列，置信水平(1-p)

殘差qq圖，QQ圖中的點的分布非常接近一條直線，只有在兩段的時候有一些偏差，大部分數據分布很好。如果這樣，該QQ圖說明殘差有很好的正態性。

qqnorm(xtimeseriesforecasts $residuals)

qqline(xtimeseriesforecasts $residuals)

# 相關圖顯示出在3階接近顯著(置信)邊界,1、4、9階超出顯著(置信)邊界，而且Ljung-Box檢驗的p值為0.01389<0.05，說明殘差沒有很好的獨立性，不符合模型參數檢驗的第二個條件。

# 為了調查預測誤差是否是平均值為零且方差為常數的正態分布(服從零均值、方差不變的正態分布)，我們可以做預測誤差的時間曲線圖和直方圖(具有正態分布曲線)：

plot.ts(xtimeseriesforecasts$residuals) # make time plot of forecast errors

plotForecastErrors(kingstimeseriesforecasts$residuals)

forecast.Arima(airarima1,h=5,level=c(99.5))

PS:參考