量化投資學習筆記11——關於時間序列你所能作的一切

關於時間序列你所能作的一切
Siddharth Yadav
翻譯自https://www.kaggle.com/thebrownviking20/everything-you-can-do-with-a-time-series
數據文件也在上面連接裏。或者上個人github代碼庫:https://github.com/zwdnet/MyQuant/tree/master/11
目標
從我註冊這個平臺的第一週，我就被時間序列分析這個主題給迷住了。本文是關於時間序列分析的許多普遍的話題的一個集合體。我寫做本文的目的是爲時間序列分析初學者和有經驗的人提供一個基本的參考。
一些重要的事情
1.本教程還在完成中，因此你每次打開它都有可能會發現有更新的內容。
2.我在寫這篇教程時已經學習過不少這個領域的課程，我還在繼續學習更多的更高級的課程以得到更多的知識和內容。
3.若是您有任何建議或者有任何主題但願本教材覆蓋，請在評論區留言。
4.若是您欣賞本文，請必定點贊（按喜歡按鈕）。這樣它能對社區有更大的意義和幫助。
首先導入相關的庫
'''python
import os
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
%matplotlib inline
from pylab import rcParams
from plotly import tools
import plotly.plotly as py
from plotly.offline import init_notebook_mode, iplot
init_notebook_mode(connected=True)
import plotly.graph_objs as go
import plotly.figure_factory as ff
import statsmodels.api as sm
from numpy.random import normal, seed
from scipy.stats import norm
from statsmodels.tsa.arima_model import ARMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.arima_model import ARIMA
import math
from sklearn.metrics import mean_squared_error
print(os.listdir("../input"))
'''

輸出
['historical-hourly-weather-data', 'stock-time-series-20050101-to-20171231']

目錄
1.介紹日期與時間
1.1 導入時間序列數據
1.2 時間序列數據的清洗與準備
1.3 數據的可視化
1.4 時間戳和週期
1.5 使用date_range
1.6 使用to_datetime
1.7 轉換與延遲(shifting and lags)
1.8 重取樣
2.金融與統計學
2.1 改變的百分率
2.2 證券收益
2.3 相繼列的絕對改變(Absolute change in sucessive rows)
2.4 比較兩個或更多的時間序列
2.5 窗口函數
2.6 OHLC圖
2.7 蠟燭圖
2.8 自相關與部分自相關
3.時間序列分解與隨機行走
3.1 趨勢、季節性和噪音
3.2 白噪音
3.3 隨機行走
3.4 穩定性(Stationarity)
4.使用statsmodels建模
4.1 AR模型
4.2 MA模型
4.3 ARMA模型
4.4 ARIMA模型
4.5 VAR模型
4.6 狀態空間模型
4.6.1 SARIMA模型
4.6.2 未觀察到的部分(Unobserved omponents)
4.6.3 動態因子模型

1.介紹日期與時間
1.1 導入時間序列數據
如何導入數據？
首先，咱們導入本教程須要的全部數據集。所需的時間序列數據的列做爲日期時間使用parse_dates參數導入，另外可使用dateframe的index_col參數來選擇索引。
咱們將使用的數據包括:
1.谷歌股票數據
2.世界各個城市的溫度數據
3.微軟股票數據
4.世界各個城市的氣壓數據
導入數據
'''python
google = pd.read_csv("input/stock-time-series-20050101-to-20171231/GOOGL_2006-01-01_to_2018-01-01.csv", index_col = "Date", parse_dates = ["Date"])
print(google.head())
'''

'''python
humidity = pd.read_csv("input/historical-hourly-weather-data/humidity.csv", index_col = "datetime", parse_dates = ['datetime'])
print(humidity.tail())
'''

1.2 時間序列數據的清洗與準備
如何準備數據?
谷歌股票數據沒有缺失項，而氣溫數據有缺失數據。使用fillna()方法，其ffill參數採用最近的有效觀測值來填充缺失值。
填充缺失值
'''python
humidity = humidity.iloc[1:]
humidity = humidity.fillna(method = "ffill")
print(humidity.head())

'''

1.3 數據集的可視化
數據集的可視化
'''python
fig = plt.figure()
humidity["Kansas City"].asfreq("M").plot()
plt.title("Humidity in Kansas City over time(Monthly frequency)")
fig.savefig("Kansas_humidity.png")
'''

'''python
fig = plt.figure()
google["2008":"2010"].plot(subplots = True, figsize = (10, 12))
plt.title("Google stocks from 2008 to 2010")
plt.savefig("google.png")
'''

1.4 時間戳和週期
什麼是時間戳和週期?它們如何有用?
時間戳是用來表明時間中的一個點。週期是時間的一段間隔。週期能夠用來檢查某一特定事件是否發生在給定的期間內。它們也能夠被轉化爲其它形式。
時間戳
'''python
timestamp = pd.Timestamp(2017, 1, 1, 12)
print(timestamp
'''

2017-01-01 12:00:00
創建一個週期
'''python

period = pd.Period("2017-01-01")
print(period)
'''

2017-01-01
檢查一個給定的時間戳是否在一個給定的時間週期中
'''python
print(period.start_time < timestamp < period.end_time)
'''

True
將時間戳轉換爲週期
'''python
new_period = timestamp.to_period(freq = "H")
print(new_period)
'''

2017-01-01 12:00
將週期轉換爲時間戳
'''python
new_timestamp = period.to_timestamp(freq = "H", how = "start")
print(new_timestamp)
'''

2017-01-01 00:00:00
1.5 使用date_range
什麼是date_range以及其爲什麼那麼有用?
date_range是一個返回固定頻率的日期時間的方法。在你基於一個已經存在的數據序列創建時間序列數據，或者從新安排整個時間序列數據時它很是有用。
以天天的頻率創建一個時間日期索引
'''python
dr1 = pd.date_range(start = "1/1/18", end = "1/9/18")
print(dr1)
'''

DatetimeIndex(['2018-01-01', '2018-01-02', '2018-01-03', '2018-01-04', '2018-01-05', '2018-01-06', '2018-01-07', '2018-01-08', '2018-01-09'], dtype='datetime64[ns]', freq='D')

以每個月的頻率創建一個時間日期索引
'''python
dr2 = pd.date_range(start = "1/1/18", end = "1/1/19", freq = "M")
print(dr2)
'''

DatetimeIndex(['2018-01-31', '2018-02-28', '2018-03-31', '2018-04-30', '2018-05-31', '2018-06-30', '2018-07-31', '2018-08-31', '2018-09-30', '2018-10-31', '2018-11-30', '2018-12-31'], dtype='datetime64[ns]', freq='M')

不設置日期起點，設定終點和週期
'''python
dr3 = pd.date_range(end = "1/4/2014", periods = 8)
print(dr3)
'''

DatetimeIndex(['2013-12-28', '2013-12-29', '2013-12-30', '2013-12-31', '2014-01-01', '2014-01-02', '2014-01-03', '2014-01-04'], dtype='datetime64[ns]', freq='D')

指定起止日期和週期
'''python

dr4 = pd.date_range(start = "2013-04-24", end = "2014-11-27", periods = 3)
print(dr4)
'''

DatetimeIndex(['2013-04-24', '2014-02-09', '2014-11-27'], dtype='datetime64[ns]', freq=None)
1.6 使用to_datetime
pandas.to_datetime()用來將變量轉換爲datetime變量。這裏，將一個DateFrame轉換爲datetime序列。
使用to_datetime
'''python
df = pd.DataFrame({
"year" : [2015, 2016],
"month" : [2, 3],
"day" : [4, 5]
})
print(df)
df = pd.to_datetime(df)
print(df)
df = pd.to_datetime("01-01-2017")
print(df)
'''

year month day
0 2015 2 4
1 2016 3 5
0 2015-02-04
1 2016-03-05
dtype: datetime64[ns]
2017-01-01 00:00:00

1.7 變換和延遲
咱們能夠經過提供時間間隔來變換索引。這對於比較一個時間序列與其自身的歷史數據頗有用。
索引變換
'''python
fig = plt.figure()
humidity["Vancouver"].asfreq('M').plot(legend = True)
shifted = humidity["Vancouver"].asfreq('M').shift(10).plot(legend = True)
shifted.legend(['Vancouver','Vancouver_lagged'])
fig.savefig("shifted.png")
'''

1.8 重取樣
上採樣(Upsampling):時間序列數據從低時間頻率到高時間頻率採樣(月到天)。這涉及到採用填充或插值的方法處理缺失數據。
下采樣(Downsampling):時間序列數據從高時間頻率到低時間頻率採樣(天到月)。這涉及到合併已存在的數據。
採用氣壓數據演示重採樣
'''python
pressure = pd.read_csv("input/historical-hourly-weather-data/pressure.csv", index_col = "datetime", parse_dates = ["datetime"])
print(pressure.tail())
pressure = pressure.iloc[1:]
'''

用前值填充nan
'''python
pressure = pressure.fillna(method = "ffill")
print(pressure.tail())
pressure = pressure.fillna(method = "bfill")
print(pressure.head())
'''

首先咱們使用ffill參數用nan以前最後一個可用數據來填充，接着咱們使用bfill用nan以後第一個可用的數據來填充。
輸出數據規模
'''python
print(pressure.shape)
'''

(45252, 36)

使用平均數從小時數據到3天數據進行向下採樣
'''python

pressure = pressure.resample("3D").mean()
print(pressure.head())
print(pressure.shape)
'''

Vancouver ... Jerusalem
datetime ...
2012-10-01 931.627119 ... 990.525424
2012-10-04 1019.083333 ... 990.083333 2012-10-07 1013.930556 ... 989.833333 2012-10-10 1015.000000 ... 987.888889 2012-10-13 1008.152778 ... 990.430556
[5 rows x 36 columns]
(629, 36)

只剩下較少的行數了。如今咱們從3天數據向每日數據進行上採樣。
從三日數據向每日數據進行上採樣
'''python
pressure = pressure.resample('D').pad()
print(pressure.head())
print(pressure.shape)
'''

Vancouver ... Jerusalem datetime ...
2012-10-01 931.627119 ... 990.525424 2012-10-02 931.627119 ... 990.525424 2012-10-03 931.627119 ... 990.525424 2012-10-04 1019.083333 ... 990.083333 2012-10-05 1019.083333 ... 990.083333
[5 rows x 36 columns]
(1885, 36)

2.金融和統計學
2.1 改變的百分率
改變的百分率
'''python
fig = plt.figure()
google["Change"] = google.High.div(google.High.shift())
google["Change"].plot(figsize = (20, 8))
fig.savefig("percent.png")
'''

2.2 證券收益
證券收益
'''python
fig = plt.figure()
google["Return"] = google.Change.sub(1).mul(100)
google["Return"].plot(figsize = (20, 8))
fig.savefig("Return1.png")
'''

另外一個計算方法
'''python
fig = plt.figure()
google.High.pct_change().mul(100).plot(figsize = (20, 6))
fig.savefig("Return2.png")
'''

2.3 相繼列的絕對改變(Absolute change in sucessive rows)
比較相繼序列的絕對改變
'''python
fig = plt.figure()
google.High.diff().plot(figsize = (20, 6))
fig.savefig("AbsoluteChange.png")
'''

2.4 比較兩個或更多時間序列
咱們將經過正態化(normalizing)來比較兩個時間序列。這是經過將時間序列中的每一個元素除以第一個元素來實現的。這樣兩個時間序列都在同一個起點開始，能夠更容易的比較。
比較兩個不一樣的序列，微軟和谷歌的股票
'''python
microsoft = pd.read_csv("input/stock-time-series-20050101-to-20171231/MSFT_2006-01-01_to_2018-01-01.csv", index_col = "Date", parse_dates = ["Date"])
'''

在正態化之前繪圖
'''python
fig = plt.figure()
google.High.plot()
microsoft.High.plot()
plt.legend(["Google", "Microsoft"])
fig.savefig("Compare.png")
'''

進行正態化並進行比較
'''python
normalized_google = google.High.div(google.High.iloc[0]).mul(100)
normalized_microsoft = microsoft.High.div(microsoft.High.iloc[0]).mul(100)
fig = plt.figure()
normalized_google.plot()
normalized_microsoft.plot()
plt.legend(["Google", "Microsoft"])
fig.savefig("NormalizedCompare.png")
'''

能夠看到谷歌的股價表現好於微軟數倍。
2.5 窗口函數
窗口函數用於定義子序列的週期，計算子週期內的子集。
有兩種:
Rolling 相同的大小和切片
Expanding 包含全部以前的數據
Rolling窗口函數
90日均線吧
'''python
rolling_google = google.High.rolling("90D").mean()
fig = plt.figure()
google.High.plot()
rolling_google.plot()
plt.legend(["High", "Rolling Mean"])
fig.savefig("RollongGoogle.png")
'''

Expanding窗口函數
'''python
microsoft_mean = microsoft.High.expanding().mean()
microsoft_std = microsoft.High.expanding().std()
fig = plt.figure()
microsoft.High.plot()
microsoft_mean.plot()
microsoft_std.plot()
plt.legend(["High", "Expanding Mean", "Expanding Standard Deviation"])
fig.savefig("ExpandingMicrosoft.png")
'''

2.6 收市價圖(OHLC charts)
收市價圖是表示必定時間週期內任何類型的價格的開盤價、最高價、最低價、以及收盤價的圖形。開盤(Open)-最高(High)-最低(Low)-收盤(Close)，OHLC圖用來做爲一種交易工具來可視化和分析證券、外匯、股票、債券、期貨等的價格隨時間的變化。收市價圖對於解釋市場價格的每日的變化以及經過模式識別來預測將來的價格改變頗有幫助。
收市價圖的y軸用來表示價格尺度，而x軸用來表示時間尺度。在每一單獨的時間週期內，一個蠟燭圖用一個符號來表明兩個範圍：交易的最高價和最低價，以及那個時間段（例如一天）的開盤價和收盤價。在符號的範圍內，最高價和最低價的範圍用主要豎線的長度來表明。開盤價和收盤價用位於豎線左邊（表明開盤價）和右邊（表明收盤價）的刻度線來表示。
每一個收市價圖符號都有顏色，以區別是「牛市」(bullish)（收盤價比開盤價高）或者「熊市」（bearish）（收盤價比開盤價低）。（文中顏色貌似與咱們的習慣是反着的，熊市是紅色，牛市是綠色。——譯者注）

(做圖這部分程序我手機上庫用不了，就照搬原文了。）
OHLC chart of June 2008
'''python
trace = go.Ohlc(x=google['06-2008'].index,
open=google['06-2008'].Open,
high=google['06-2008'].High,
low=google['06-2008'].Low,
close=google['06-2008'].Close)
data = [trace]iplot(data, filename='simple_ohlc')
'''

OHLC chart of 2008
'''python
trace = go.Ohlc(x=google['2008'].index,
open=google['2008'].Open,
high=google['2008'].High,
low=google['2008'].Low,
close=google['2008'].Close)
data = [trace]iplot(data, filename='simple_ohlc')
'''

OHLC chart of 2008
'''python
trace = go.Ohlc(x=google.index,
open=google.Open,
high=google.High,
low=google.Low,
close=google.Close)
data = [trace]iplot(data, filename='simple_ohlc')
'''

2.7 蠟燭圖
這種圖用來做爲一種交易工具可視化分析證券、衍生品、外匯、股票、債券、期貨等的價格隨時間的運動。儘管蠟燭圖使用的符號像一個箱子的符號，可是它們的功能是不一樣的，不能彼此混淆。
蠟燭圖使用一個蠟燭樣的符號來顯示多種價格信息，例如開盤價、收盤價、最高價和最低價等。每一個符號表明一個單獨時間間隔內的壓縮的交易活動（一分鐘、一小時、一天、一個月等）。每一個蠟燭符號畫在x軸上一個單獨的時間尺度上，以顯示那段時間的交易活動。
符號中間的矩形被稱爲實體，用來顯示那段時間的開盤價與收盤價的範圍。從其頂部和底部延伸出來的線段表明下影和上影（或者燈芯）。它們表明那段時間裏的最低價和最高價。當爲牛市時（收盤價高於開盤價），實體常爲白色或綠色。當爲熊市時（收盤價低於開盤價），實體被塗爲黑色或紅色。(仍是跟咱們的習慣相反）

經過其顏色和形狀，蠟燭圖對於探測和預測市場隨時間的趨勢以及解釋市場的每日的變更是很是好的工具。例如，實體越長，賣壓或買壓約大。實體越短，意味着該時間短內價格變更很是小。
蠟燭圖經過各類指標幫助顯示市場心理（交易者的恐懼或貪婪），如形狀和顏色，還有一些能夠從蠟燭圖中識別出來的模式。總的來講，大約有42種已經被識別出來的或簡單或複雜的模式。這些從蠟燭圖中識別出來的模式對於顯示價格關係和預測市場將來可能的變更頗有幫助。這裏有一些模式。

請注意，蠟燭圖並不表示在開盤價格和收盤價格之間發生的事件——只是關於兩種價格之間的關係。所以你沒法指出這個時間斷內的交易的波動性。
來源： https://datavizcatalogue.com/methods/candlestick_chart.html
Candlestick chart of march 2008
'''python
trace = go.Candlestick(x=google['03-2008'].index,
open=google['03-2008'].Open,
high=google['03-2008'].High,
low=google['03-2008'].Low,
close=google['03-2008'].Close)
data = [trace]iplot(data, filename='simple_candlestick')
'''

Candlestick chart of 2008
'''python
trace = go.Candlestick(x=google['2008'].index,
open=google['2008'].Open,
high=google['2008'].High,
low=google['2008'].Low,
close=google['2008'].Close)
data = [trace]iplot(data, filename='simple_candlestick')
'''

Candlestick chart of 2006-2018
'''python
trace = go.Candlestick(x=google.index,
open=google.Open,
high=google.High,
low=google.Low,
close=google.Close)
data = [trace]iplot(data, filename='simple_candlestick')
'''

2.8 自相關與部分自相關
自相關(Autocorrelation)——自相關函數(The autocorrelation function, ACF)測量一個序列在不一樣的片斷上與其自身的相關性。
部分自相關(Partial Autocorrelation)——部分自相關函數能夠被解釋爲一個序列的片斷是其以前的片斷的迴歸。能夠用標準線性迴歸來解釋這個概念，這是部分片斷的改變而其它片斷保持不變。
來源：https://www.quora.com/What-is-the-difference-among-auto-correlation-partial-auto-correlation-and-inverse-auto-correlation-while-modelling-an-ARIMA-series
自相關性
自相關
'''python
fig = plt.figure()
plot_acf(humidity["San Diego"], lags = 25, title = "San Diego")
plt.savefig("acf.png")
'''

全部的延遲都接近1或者至少大於置信區間，它們有顯著的統計學差別的。
部分自相關
部分自相關
'''python
fig = plt.figure()
plot_pacf(humidity["San Diego"], lags = 25, title = "San Diego, pacf")
plt.savefig("pacf.png")
'''

因爲差別有統計學意義，在頭兩個延遲以後的部分自相關性很是低。
'''python
plot_pacf(microsoft["Close"], lags = 25)
plt.savefig("ms_pacf.png")
'''

這裏，只有第0,1和20個延遲差別有統計學意義。

3.時間序列分解與隨機行走
3.1 趨勢，季節性和噪音
這些是一個時間序列的組成部分
趨勢：一個時間序列的持續向上或向下的斜率
季節性：一個時間序列的清晰的週期性模式（就像正弦函數）
噪音：異常或缺失數據
趨勢，季節性和噪音
'''python
fig = plt.figure()
google["High"].plot(figsize = (16, 8))
fig.savefig("google_trend.png")
'''

分解
'''python
decomposed_google_volume = sm.tsa.seasonal_decompose(google["High"], freq = 360)
fig = decomposed_google_volume.plot()
fig.savefig("decomposed.png")
'''

在上圖中有一個清晰的向上的趨勢。
你也能夠看到規則的週期性改變
不規則的噪音表明着數據異常或缺失。

3.2 白噪音
白噪音是
恆定的平均值
恆定的差別
全部的偏移的零自相關
白噪音
'''python
fig = plt.figure()
rcParams["figure.figsize"] = 16, 6
white_noise = np.random.normal(loc = 0, scale = 1, size = 1000)
plt.plot(white_noise)
fig.savefig("whitenoise.png")
'''

繪製白噪音的自相關關係

'''python
fig = plt.figure()
plot_acf(white_noise, lags = 20)
plt.savefig("wn_acf.png")
'''

觀察全部偏移都在置信區間內（陰影部分），所以是有統計學意義的。

3.3 隨機行走
隨機行走是一個數學概念，是一個隨機的過程，描述一個在一些數學空間（如整數）由持續的隨機步數來描述的路徑。
通常來講，對於股市，今天的價格 = 昨天的價格+噪音
Pt = Pt-1 + εt
隨機行走不能被預測，由於噪音是隨機的。
具備Drift(drift(μ) 是 0-平均值）的隨機行走
Pt - Pt-1 = μ + εt
對隨機行走的迴歸測試
Pt = α + βPt-1 + εt
化簡爲
Pt - Pt-1 = α + βPt-1 + εt
檢驗：
H0：β = 1（這是一個隨機行走）
H1：β < 1 （這不是一個隨機行走）
Dickey-Fuller（DF）檢驗
H0：β = 0（這是一個隨機行走）
H1：β < 0 （這不是一個隨機行走）
單位根檢驗（Augmented Dickey-Fuller test）
單位根檢驗的零假設是一個時間序列樣本存在單位根，基本上單位根檢驗在RHS上有更多的延遲改變。

單位根檢驗谷歌和微軟的成交量

'''python
adf = adfuller(microsoft["Volume"])
print("p-value of microsoft: {}".format(float(adf[1])))
adf = adfuller(google["Volume"])
print("p-value of google: {}".format(float(adf[1])))
'''

p-value of microsoft: 0.0003201525277652296
p-value of google: 6.510719605767195e-07
因爲微軟的p值爲0.0003201525小於0.05，零假設被拒絕，這個序列不是一個隨機行走序列。
谷歌的p值爲0.0000006510小於0.05（此處原文爲大於，疑有誤），零假設被拒絕，這個序列不是一個隨機行走序列。
產生一個隨機行走序列
'''python
fig = plt.figure()
seed(42)
rcParams["figure.figsize"] = 16, 6
random_walk = normal(loc = 0, scale = 0.01, size = 1000)
plt.plot(random_walk)
fig.savefig("random_walk.png")
'''

'''python
fig = plt.figure()
plt.hist(random_walk)
fig.savefig("random_hist.png")
'''

3.4 穩定性
一個穩定的時間序列是指其統計性質如平均值，方差，自相關性等，都不隨時間變化的時間序列。
強穩定性：是指其機率分佈絕對的不隨時間變化而變化的隨機過程。所以，相似平均值、方差等參數也不隨時間而變化。
弱穩定性：是指平均值、方差、自相關性都隨時間變化保持恆定的過程。
穩定性很是重要，由於非穩定序列有不少影響因素，在建模的時候會很複雜。diff()方法能夠容易的將一個非穩定序列轉化爲穩定序列。
咱們將嘗試將上面的時間序列的週期性部分分解。（We will try to decompose seasonal component of the above decomposed time series.）
初始的非穩定序列
'''python
fig = plt.figure()
decomposed_google_volume.trend.plot()
fig.savefig("nonstationary.png")
'''

新的穩定的序列，即一階差分

'''python
fig = plt.figure()
decomposed_google_volume.trend.diff().plot()
fig.savefig("stationary.png")
'''

4.使用statstools建模
4.1 AR模型
一個自迴歸(autoregressive, AR）模型表明了一類隨機過程，這類過程用於描述天然界中特定的時間變化過程，如經濟學。該模型認爲輸出變量線性的依賴於其以前的數據以及一個隨機成分（一個有缺陷的預測值）；所以模型是以隨機差分方程的形式出現。
AR(1)模型
Rt = μ + ΦRt-1 + εt
因爲RHS只有一個延遲值(Rt-1)，這被稱爲一階AR模型，其中μ是平均值， ε是t時刻的噪音。
若是Φ=1，這是隨機行走。若是Φ=0，這是噪音。若是-1<Φ<1，它是穩定的。若是Φ爲負值，有一我的爲因素，若是Φ爲正值，有一個動量。( If ϕ is -ve, there is men reversion. If ϕ is +ve, there is momentum.)
AR(2)模型
Rt = μ + Φ1Rt-1 + Φ2Rt-2 + εt
AR(3)模型
Rt = μ + Φ1Rt-1 + Φ2Rt-2 + Φ3Rt-3 + εt
AR(1)模擬模型
AR(1) MA(1)模型: AR參數 = 0.9

'''python
fig = plt.figure()
rcParams['figure.figsize'] = 16, 12
plt.subplot(4,1,1)
ar1 = np.array([1, -0.9])
ma1 = np.array([1])
AR1 = ArmaProcess(ar1, ma1)
sim1 = AR1.generate_sample(nsample = 1000)
plt.title("AR(1) model : AR parameter = +0.9")
plt.plot(sim1)
'''

AR(1) MA(1)模型: AR參數 = -0.9
'''python
plt.subplot(4,1,2)
ar2 = np.array([1, 0.9])
ma2 = np.array([1])
AR2 = ArmaProcess(ar2, ma2)
sim2 = AR2.generate_sample(nsample = 1000)
plt.title("AR(1) model : AR parameter = -0.9")
plt.plot(sim2)
'''

AR(2) MA(1)模型: AR參數 = 0.9
'''python
plt.subplot(4,1,3)
ar3 = np.array([2, -0.9])
ma3 = np.array([1])
AR3 = ArmaProcess(ar3, ma3)
sim3 = AR3.generate_sample(nsample = 1000)
plt.title("AR(2) model : AR parameter = +0.9")
plt.plot(sim3)
'''

AR(2) MA(1)模型: AR參數 = -0.9
'''python
plt.subplot(4,1,4)
ar4 = np.array([2, 0.9])
ma4 = np.array([1])
AR4 = ArmaProcess(ar4, ma4)
sim4 = AR4.generate_sample(nsample = 1000)
plt.title("AR(2) model : AR parameter = -0.9")
plt.plot(sim4)
fig.savefig("AR.png")
'''

一個模擬模型的預測
'''python
model = ARMA(sim1, order=(1, 0))
result = model.fit()
print(result.summary())
print("μ = {}, φ = {}".format(result.params[0], result.params[1]))
'''

Φ約爲0.9，是咱們在第一個模擬模型中選擇的AR參數。
預測模型
'''python
fig = plt.figure()
fig = result.plot_predict(start = 900, end = 1010)
fig.savefig("AR_predict.png")
'''

'''python
rmse = math.sqrt(mean_squared_error(sim1[900:1011], result.predict(start = 900, end = 999)))
print("The root mean squared error is {}.".format(rmse))
'''

The root mean squared error is 1.0408054544358292.
y的預測值已經畫出，很整潔！
預測蒙特利爾的溼度
'''python

humid = ARMA(humidity["Montreal"].diff().iloc[1:].values, order = (1, 0))
res = humid.fit()
fig = plt.figure()
fig = res.plot_predict(start = 1000, end = 1100)
fig.savefig("humid_arma.png")
'''

'''python

rmse = math.sqrt(mean_squared_error(humidity["Montreal"].diff().iloc[900:1000].values, result.predict(start=900,end=999)))print("The root mean squared error is {}.".format(rmse))
'''

The root mean squared error is 7.218388589479766.
不是很使人印象深入，但讓咱們試試谷歌股票。
預測谷歌的收盤價
'''python

humid = ARMA(google["Close"].diff().iloc[1:].values, order = (1, 0))
res = humid.fit()
fig = plt.figure()
fig = res.plot_predict(start = 900, end = 1100)
fig.savefig("google_arma.png")
'''

總有更好的模型

4.2 MA模型
移動平均模型在單變量時間序列中頗有用。模型假設輸出變量線性的依賴於當前的各類隨後的隨機變量（不許確預測值）。
MA(1)模型
Rt = μ + εt1 + θεt-1
即今日的收益= 平均值+今日的噪音+昨日的噪音。
由於RHS中只有一個延遲，這是一階的MA模型。
MA(1)模擬模型
rcParams["figure.figsize"] = 16, 6
ar1 = np.array([1])
ma1 = np.array([1, -0.5])
MA1 = ArmaProcess(ar1, ma1)
sim1 = MA1.generate_sample(nsample = 1000)
plt.plot(sim1)
plt.savefig("ma1.png")

MA模型的建模
創建MA模型的預測
model = ARMA(sim1, order=(0, 1))
result = model.fit()
print(result.summary())
print("μ={} ,θ={}".format(result.params[0],result.params[1]))

μ=-0.02284716848276931 ,θ=-0.5650012559991154
MA模型的預測
使用MA模型進行預測
model = ARMA(humidity["Montreal"].diff().iloc[1:].values, order=(0, 3))
result = model.fit()
print(result.summary())
print("μ={} ,θ={}".format(result.params[0],result.params[1]))
result.plot_predict(start = 1000, end = 1100)
plt.savefig("ma_forcast.png")

rmse = math.sqrt(mean_squared_error(humidity["Montreal"].diff().iloc[1000:1101].values, result.predict(start=1000,end=1100)))
print("The root mean squared error is {}.".format(rmse))
The root mean squared error is 11.345129665763626.
接着，是ARMA模型
4.3 ARMA模型
自迴歸移動平均模型（Autoregressive–moving-average，ARMA）提供一個以二項式形式描述一個（弱的）穩定隨機過程的模型。一個是自迴歸，另外一個是移動平均。它是AR和MA模型的綜合。
ARMA(1,1)模型
Rt = μ + ΦRt-1 + εt + θεt-1
基本上，它表明着今日收益 = 平均值 + 昨日的收益 + 噪音 + 昨日的噪音。
ARMA預測模型的建模
由於與AR和MA模型相似，就不進行模擬了。直接進行預測。
模擬和預測微軟股票的市值
model = ARMA(microsoft["Volume"].diff().iloc[1:].values, order = (3, 3))
result = model.fit()
print(result.summary())
print("μ={} ,θ={}".format(result.params[0],result.params[1]))
result.plot_predict(start = 1000, end = 1100)
plt.savefig("arma_forcast.png")

rmse = math.sqrt(mean_squared_error(microsoft["Volume"].diff().iloc[1000:1101].values, result.predict(start=1000,end=1100)))
print("The root mean squared error is {}.".format(rmse))
The root mean squared error is 38038241.66905847.
ARMA模型的預測結果要優於AR和MA模型

4.4 ARIMA模型
求和自迴歸移動平均模型（autoregressive integrated moving average ，ARIMA）是ARMA模型的通常化。這些模型都是都是擬合時間序列數據，以便更好的理解數據或者預測將來的數據。它應用在不穩定的序列數據，經過一系列的差分步驟（模型相應的「求和」部分）能夠消除數據的不穩定。ARIMA模型以ARIMA(p, d, q)形式表示，p是AR的參數，d是差分參數，q是MA的參數
ARIMA(1, 0, 0)
yt = a1yt-1 + εt

ARIMA(1, 0, 1)
yt = a1yt-1 + εt + b1εt-1

ARIMA(1, 1, 1)
Δyt = a1Δyt-1 + εt + b1εt-1, 其中Δyt = yt - yt-1

創建ARIMA的預測模型
使用ARIMA模型進行預測
預測微軟股票的市值
rcParams["figure.figsize"] = 16, 6
model = ARIMA(microsoft["Volume"].diff().iloc[1:].values, order = (2, 1, 0))
result = model.fit()
print(result.summary())
result.plot_predict(start = 700, end = 1000)
plt.savefig("Arima_predict.png")

rmse = math.sqrt(mean_squared_error(microsoft["Volume"].diff().iloc[700:1001].values, result.predict(start=700,end=1000)))
print("The root mean squared error is {}.".format(rmse))

The root mean squared error is 61937593.98493614.
考慮一個輕微的延遲，這是一個很好的模型。
4.5 VAR模型
向量自迴歸(Vector autoregression, VAR)是一個隨機過程模型，用來捕捉多個時間序列之間的線性相關性。VAR模型是單變量自迴歸模型(AR模型)推廣到多個變量的狀況。在VAR模型中全部變量進入模型的途徑都一致：每一個變量都有一個方程基於其本身的延遲值，其它模型變量的延遲值，以及一個偏差因子來解釋其演變。VAR模型不須要更多的關於影響一個變量的因素的知識，就像在結構化模型中那樣：模型須要的惟一的先導知識是變量列表，其中的變量被暫時地假設會彼此相互影響。
VAR模型
預測谷歌和微軟的收盤價
train_sample = pd.concat([google["Close"].diff().iloc[1:], microsoft["Close"].diff().iloc[1:]], axis=1)
model = sm.tsa.VARMAX(train_sample, order = (2, 1), trend = 'c')
result = model.fit(maxiter = 1000, disp = False)
print(result.summary())
predicted_result = result.predict(start = 0, end = 1000)
fig = result.plot_diagnostics()
fig.savefig("Var_predict.png")
計算偏差
rmse = math.sqrt(mean_squared_error(train_sample.iloc[1:1002].values, predicted_result.values))
print("The root mean squared error is {}.".format(rmse))

4.6 狀態空間模型
一個通常的狀態空間模型的形式爲：
yt = Ztαt + dt + εt
αt = Ttαt - 1 + ct + Rtηt
其中yt表明了在時間t中的觀察向量，αt表明了（未觀察的）在時間t的狀態向量，不規則的成分定義以下：
εt ~ N(0, Ht)
ηt ~ N(0, Qt)
方程中的其他變量(Zt, dt, Ht, Tt, ct, Rt, Qt)是描述過程的矩陣。它們的名稱和維度以下：
(這些就不翻譯了)
Z : design (k_endog×k_states×nobs)
d : obs_intercept (k_endog×nobs)
H : obs_cov (k_endog×k_endog×nobs)
T : transition (k_states×k_states×nobs)
c : state_intercept (k_states×nobs)
R : selection (k_states×k_posdef×nobs)
Q : state_cov (k_posdef×k_posdef×nobs)
若是一個矩陣是時間不變的（例如，Zt = Zt + 1 ∀ t），其最後的維度也許大小爲1而不是節點的數量。
這個通常形式歸納了許多很是流行的線性時間序列模型（以下）而且有很高的擴展性，容許對缺失的觀察進行估計，進行預測，推進響應函數，等等。
來源:https://www.statsmodels.org/dev/statespace.html
4.6.1 SARIMA模型
SARIMA模型對於季節性時間序列的建模頗有用，其數據的平均值和其它統計指標在某一年度內是不穩定的，SARIMA模型是非季節性自迴歸移動平均模型(ARMA)和自迴歸求和移動平均模型(ARIMA）的直接擴展。
SARIMA模型
預測谷歌的收盤價
train_sample = google["Close"].diff().iloc[1:].values
model = sm.tsa.SARIMAX(train_sample, order = (4, 0, 4), trend = 'c')
result = model.fit(maxiter = 1000, disp = True)
print(result.summary())
predicted_result = result.predict(start = 0, end = 500)
fig = result.plot_diagnostics()
fig.savefig("sarimax.png")
計算偏差
rmse = math.sqrt(mean_squared_error(train_sample[1:502], predicted_result))
print("The root mean squared error is {}.".format(rmse))

fig = plt.figure()
plt.plot(train_sample[1:502], color = "red")
plt.plot(predicted_result, color = "blue")
plt.legend(["Actual", "Predicted"])
plt.title("Google closing price")
fig.savefig("sarimax_test.png")

4.6.2 未觀察成分
UCM將序列拆分爲組成成分，例如趨勢、季節、週期，以及衰退因素，以預測序列。下面的模型給出了一個可能的方案。

來源：http://support.sas.com/documentation/cdl/en/etsug/66840/HTML/default/viewer.htm#etsug_ucm_details01.htm
未觀察成分模型
預測谷歌的收盤價
train_sample = google["Close"].diff().iloc[1:].values
model = sm.tsa.UnobservedComponents(train_sample, "local level")
result = model.fit(maxiter = 1000, disp = True)
print(result.summary())
predicted_result = result.predict(start = 0, end = 500)
fig = result.plot_diagnostics()
fig.savefig("unobserve.png")
計算偏差
rmse = math.sqrt(mean_squared_error(train_sample[1:502], predicted_result))
print("The root mean squared error is {}.".format(rmse))
fig = plt.figure()
plt.plot(train_sample[1:502], color = "red")
plt.plot(predicted_result, color = "blue")
plt.legend(["Actual", "Predicted"])
plt.title("Google closing price")
fig.savefig("unobserve_test.png")

4.6.3 動態因子模型
動態因子模型是一個靈活的模型，用於多變量時間序列，其中觀測的內部變量與外部協變量和未觀測的因子呈線性關係。所以有一個向量自迴歸結構。未觀測的因子也許是外部協變量的一個函數。依賴變量對方程的干擾也許是自相關的。
動態因子模型
預測谷歌的收盤價
train_sample = pd.concat([google["Close"].diff().iloc[1:], microsoft["Close"].diff().iloc[1:]], axis=1)
model = sm.tsa.DynamicFactor(train_sample, k_factors=1, factor_order=2)
result = model.fit(maxiter = 1000, disp = True)
print(result.summary())
predicted_result = result.predict(start = 0, end = 1000)
fig = result.plot_diagnostics()
fig.savefig("DynamicFactor.png")
計算偏差
rmse = math.sqrt(mean_squared_error(train_sample[1:502], predicted_result))
print("The root mean squared error is {}.".format(rmse))

我將盡快增長更多的迴歸模型，覆蓋更多的主題。但根據個人經驗，對於時間序列預測的最好模型是LSTM，它是基於循環神經網絡( Recurrent Neural Networks）的。我已經爲這個主題準備了一個教程，這是連接： https://www.kaggle.com/thebrownviking20/intro-to-recurrent-neural-networks-lstm-gru

參考文獻（有更深刻的內容和解釋）：

• Manipulating Time Series Data in Python https://www.datacamp.com/courses/manipulating-time-series-data-in-python
• Introduction to Time Series Analysis in Python https://www.datacamp.com/courses/introduction-to-time-series-analysis-in-python
• Visualizing Time Series Data in Python https://www.datacamp.com/courses/visualizing-time-series-data-in-python
• VAR models and LSTM https://www.youtube.com/watch?v=_vQ0W_qXMxk
• State space models https://www.statsmodels.org/dev/statespace.html
  敬請期待更多的內容！而且別忘了點贊和評論。

(原文全文完)
本文寫於春節前，本想收假上班就發。因爲衆所周知的緣由，假期一延再延，我至今沒有肯定準確的上班時間。所以仍是發了吧。因爲是從手機上發的，格式可能會亂，尤爲代碼部分。等收假能用電腦再改了。祝疫情早日結束。祝你們身體健康。

我發文章的四個地方，歡迎你們在朋友圈等地方分享，歡迎點「在看」。
個人我的博客地址：https://zwdnet.github.io
個人知乎文章地址： https://www.zhihu.com/people/zhao-you-min/posts
個人博客園博客地址： https://www.cnblogs.com/zwdnet/ 個人微信我的訂閱號：趙瑜敏的口腔醫學學習園地