python与时间序列基本教程4(超过1.9万字 代码超过900行 包括49个图)

介绍

代码、数据全部免费，都放在我的gitee仓库里面https://gitee.com/yuanzhoulvpi/time_series，想要使用的可以直接到我这个仓库下载。

本文是继python与时间序列(开篇)的第4篇，详细介绍了一些时间序列的描述性统计知识，并且附有代码、数据。百分百可以复现。

本文内容比较多【初略统计：所有文本超过1.9万字，代码有900行】。

jupyter notebook目录如下：

• 介绍时间序列

  • 导入时间序列数据
  • 处理(清洗)数据
  • 对数据进行可视化
  • 时间戳和时间段
  • 使用date_range
  • 滞后
  • 重采样
  • 金融统计数据
  • 还手、回报、差分
  • 多个时间序列的比较
  • 窗口函数
  • OHLC曲线、烛箱图
  • 自相关和偏自相关

• 时间序列的分解

  • 白噪声
  • 随机游走

• 使用python建模

  • AR模型
  • MA模型
  • ARMA、ARIMA、SARIMA模型
  • Var模型
  • 结构时间序列模型
  • 动态因子模型

介绍时间序列

导入需要的包

# Importing libraries
import os
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
# # Above is a special style template for matplotlib, highly useful for visualizing time series data
# %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
import seaborn as sns
from sklearn.metrics import mean_squared_error

导入时间序列数据

# 在使用pandas读取文件的时候，可以使用index_col将数据的某一列直接设置为列索引；
# parse_date这个参数可以将某一列数据自动解析为时间格式。

google = pd.read_csv("../datasets/GOOGL_2006-01-01_to_2018-01-01.csv", index_col='Date', parse_dates=['Date'])
google.head(5)

#%%

humidity = pd.read_csv("../datasets/humidity.csv", index_col='datetime', parse_dates=['datetime'])
humidity.head(4)

处理(清洗)数据

上面的数据中, google数据是没有缺失值的，但是humidity数据是有缺失值的，因此，需要对数据做一下处理。

# 观察到humidity数据的第一行是空的，而且每一个时间基本上都非常接近的，因此我们将第一行给丢掉，然后使用向前填充数据

humidity = humidity.iloc[1:]
humidity = humidity.fillna(method='ffill')
humidity.head(4)

对数据进行可视化

# 将时间序列转换成特定的频率，下面的.asfreq('M')中，就是将时间序列数据转换成频率为月的
fig, ax = plt.subplots(figsize=(10,4), dpi=300)
humidity["Kansas City"].asfreq('M').plot(ax=ax, c='blue') # asfreq method is used to convert a time series to a specified frequency. Here it is monthly frequency.
plt.title('Humidity in Kansas City over time(Monthly frequency)')
plt.show()

bigcolor = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
temp_colname = ['Open','High', 'Low', 'Close', 'Volume']
fig, ax = plt.subplots(nrows=len(temp_colname), ncols=1, figsize=(10,12),
                       sharex=True,
                       dpi=300)

for index, colname_ in enumerate(temp_colname):
    tempdata = google['2008':'2010'][colname_]
    ax[index].plot(tempdata.index,tempdata, label=colname_, c=bigcolor[index])
    ax[index].legend()
    # google['2008':'2010'][colname_].plot(ax=ax[index], label=colname_)
fig.suptitle("Google stock attributes from 2008 to 2010",y=0.92)
plt.xticks(rotation=45)

时间戳和时间段

经常处理数据的可以遇到这样的单词：Timestamps，其实可以叫时间戳，本质上就是一个时间点。而Period,可以叫时间段，强调的是一段时间。

#%%

# 创建一个时间点
timestamp = pd.Timestamp(2021,10,1,12)
timestamp

#%%

# 创建一个时间段
period = pd.Period('2021-10-01')
# 这个时间段，默认指的就是这一天
period

#%%

# 看看上面的时间点是不是在这个时间段里面
period.start_time < timestamp < period.end_time

#%%

# 将时间点转换为时间段
new_period = timestamp.to_period(freq='H')
new_period

使用date_range

date_range可以返回一个固定频率的日期时间格式的索引，对已有的数据创建时间索引，然后处理数据。

#%%

# 创建一个频率以天的时间索引
dr1 = pd.date_range(start='2021-10-01' ,end='2021-10-07')
dr1


#%%

# 创建一个频率为月的时间索引
dr2 = pd.date_range(start='2020-10-01', end='2021-11-01', freq='M')
dr2

#%%

# 创建一个以特定时间周期的时间索引（这里将一个时间段分成了5个时间点）
dr3 = pd.date_range(start='2021-10-01', end='2021-11-14', periods=5)
dr3

#%%

# 如果只是同时设置start和period，会将时间向后衍生period天；
dr4 = pd.date_range(start='2021-10-01',periods=5)
dr4

#%%

# 如果只是同时设置end和period，会将时间向前衍生period天；
dr5 = pd.date_range(end='2021-11-14', periods=5)
dr5

滞后

使用pandas的shift就行

#%%

fig, ax = plt.subplots(figsize=(10,4), dpi=300)
humidity['Vancouver'].asfreq('M').plot(ax=ax, legend=True)
shifted = humidity['Vancouver'].asfreq('M').shift(10).plot(ax=ax, legend=True)
shifted.legend(['Vancouver', 'Vancouver_lagged'])
fig.suptitle(t="Vancouver_lagged_10 VS Vancouver")

重采样

1. 上采样：指的就是将数据从低频数据上升到高频数据，比如从月度数据上升到日数据。涉及到插入或者填充数据。
2. 下采样：指的就是将时间从高平数据下降到低频数据，比如从日数据下降到月数据。涉及到的就是一些聚合操作。

#%%

# 使用pressure数据来做这样的操作
pressure = pd.read_csv('../datasets/pressure.csv', index_col='datetime', parse_dates=['datetime'])
pressure.tail()


#%%

# 填充缺失数据，使用向前填充
pressure = pressure.iloc[1:]
pressure = pressure.fillna(method='ffill')
pressure.tail()


#%%

# 填充缺失数据，使用向后填充
pressure = pressure.fillna(method='bfill')
pressure.head()


#%%

# 在采样之前，查看pressure数据的大小
pressure.shape

#%%

# 开始降低采样，3天一次
pressure = pressure.resample('3D').mean()
pressure.head()

#%%

# 查看这个降采样后的数据大小
pressure.shape

#%%

# 这个时候再进行上采样
pressure = pressure.resample('D').pad()
pressure.head()


#%%

# 上采样之后，查看数据大小
pressure.shape

金融统计

查看换手率

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
google['Change'] = google.High.div(google.High.shift())
google['Change'].plot(ax=ax, linewidth=1)

回报

#%%

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
google['Return'] = google.Change.sub(1).mul(100)
google['Return'].plot(ax=ax, linewidth=1)

# 另外一种方式计算回报

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
google.High.pct_change().mul(100).plot(ax=ax, linewidth=1)

差分

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
google.High.diff().plot(ax=ax, linewidth=1)

比较两个或多个时间序列数据

# 这里选择微软的数据和谷歌数据进行比较
microsoft = pd.read_csv('../datasets/MSFT_2006-01-01_to_2018-01-01.csv', index_col='Date', parse_dates=['Date'])

#%%

# 在还没对数据进行标准化之前，看看原始数据是什么样子的

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
google.High.plot(ax=ax)
microsoft.High.plot(ax=ax)
plt.legend(['Google','Microsoft'])

# 对两个数据进行标准化，然后再进行比较
# 将数据全部除以第一个时间点的数据，然后再乘以100
fig, ax = plt.subplots(figsize=(10, 4), dpi=300)

normalized_google = google.High.div(google.High.iloc[0]).mul(100)
normalized_microsoft = microsoft.High.div(microsoft.High.iloc[0]).mul(100)
normalized_google.plot(ax=ax)
normalized_microsoft.plot(ax=ax)
plt.legend(['Google','Microsoft'])

基本上可以看出来，谷歌的增长要远远比微软增长的要快一点

窗口函数

窗口函数可以识别子周期，可以计算子周期的一些性能数据

1. Rolling:窗口大小保持不变，窗口在数据上滑动，
2. Expanding: 窗口大小不断拉长，计算的数据越来越多。

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)

# Rolling window functions
rolling_google = google.High.rolling('90D').mean()
google.High.plot(ax=ax)
rolling_google.plot(ax=ax)
plt.legend(['High','Rolling Mean'])

可以看出来，均值滚动的数据比原来的数据要更加平滑

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)

# Expanding window functions
microsoft_mean = microsoft.High.expanding().mean()
microsoft_std = microsoft.High.expanding().std()
microsoft.High.plot(ax=ax)
microsoft_mean.plot(ax=ax)
microsoft_std.plot(ax=ax)
plt.legend(['High','Expanding Mean','Expanding Standard Deviation'])

OHLC曲线

这个曲线说白了，就是显示了股票数据的OPEN、HIGH、LOW、CLOSE的数据。如果一天的OPEN是大于CLOSE那么就是红色的（因为当天跌了）；如果一天的OPEN是小于CLOSE就是绿色（因为当天涨了）

# OHLC chart of June 2008
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')
fig = go.Figure()
fig.add_trace(trace)
fig.update_layout(template='plotly_white',title='sample OHLC')

# OHLC chart of 2008
trace = go.Ohlc(x=google['2008'].index,
                open=google['2008'].Open,
                high=google['2008'].High,
                low=google['2008'].Low,
                close=google['2008'].Close)
fig = go.Figure()
fig.add_trace(trace)
fig.update_layout(template='plotly_white',title='sample OHLC')

#%%

# OHLC chart of 2008
trace = go.Ohlc(x=google.index,
                open=google.Open,
                high=google.High,
                low=google.Low,
                close=google.Close)
fig = go.Figure()
fig.add_trace(trace)
fig.update_layout(template='plotly_white',title='sample OHLC')

烛箱图

# Candlestick chart of march 2008
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)
fig = go.Figure()

fig.add_trace(trace)
fig.update_layout(template='plotly_white',title='simple_candlestick')

自相关和偏自相关

#  humidity of San Diego的自相关
fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
plot_acf(humidity["San Diego"],lags=25,title="San Diego", ax=ax)
plt.show()

由于所有滞后都接近 1 或至少大于置信区间，因此它们具有统计显着性

# humidity of San Diego的偏自相关
fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
plot_pacf(humidity["San Diego"],lags=25,ax=ax)
plt.show()

在两个滞后项以后，系数非常低了

时间序列的分解

任何时间序列可以可以被拆分为3个部分：

1. 趋势：趋势是比较明显的，比如极速的上升或者迅速下跌。
2. 季节性：可以在数据中看到明显的周期性，并且这个周期性和时间周期有关。这个周期可能是月，可能是季度，也可能是年。
3. 误差项。

但是不是说所有的时间序列都必须具备这3个部分。时间序列可能没有明显的趋势、可能没有明显的周期性。或者两个都没有。

因此，可以将时间序列看成趋势、周期性、误差项的组合。

decomposed_google_volume = sm.tsa.seasonal_decompose(google["High"],freq=360) # The frequncy is annual

# 画图
fig, ax = plt.subplots(ncols=1, nrows=4, figsize=(10, 12), sharex=True)
font = {'family': 'serif',
        'color': 'darkred',
        'weight': 'normal',
        'size': 16,
        }

def plot_decompose(result, ax,  title, fontdict=font):
    ax[0].set_title(title, fontdict=fontdict)
    result.observed.plot(ax=ax[0])
    ax[0].set_ylabel("Observed")

    result.trend.plot(ax=ax[1])
    ax[1].set_ylabel("Trend")

    result.seasonal.plot(ax=ax[2])
    ax[2].set_ylabel("Seasonal")

    result.resid.plot(ax=ax[3])
    ax[3].set_ylabel("Resid")

plot_decompose(result=decomposed_google_volume, ax=ax, title="Google High decompose")
plt.show()

可以看出来，谷歌的High有很强的趋势性、季节性、而且上面的残差还有一些模式信息没有被提取完全

白噪声

1. 均值为固定；
2. 方差固定；
3. 滞后多少项，自相关都是0；

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
white_noise = np.random.normal(loc=0, scale=1, size=1000)
ax.plot(white_noise)
plt.show()

# 查看白噪声的acf

fig, ax = plt.subplots(ncols=1,nrows=1,figsize=(10, 4), dpi=300)
plot_acf(white_noise,lags=20, ax=ax, title="white_noise acf")
plt.show()

随机游走

使用adf可以检验一个时间序列是否为随机游走

# Augmented Dickey-Fuller test on volume of google and microsoft stocks
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值都是小于0.05，拒绝原假设，所以这两个时间序列都不是随机游走。

生成一个随机游走数据

fig, ax = plt.subplots(figsize=(10, 4), dpi=300)

random_walk = normal(loc=0, scale=0.01, size=1000)
ax.plot(random_walk)
plt.show()

# 查看这个随机游走模型的分布
fig, ax = plt.subplots(figsize=(10,4), dpi=300)
sns.distplot(random_walk, hist=True, kde=True,
             bins=40, color = 'darkblue',
             hist_kws={'edgecolor':'black'},
             kde_kws={'linewidth': 4}, ax=ax)
ax.set_title("Random Walk")

使用统计工具箱

AR模型

自回归模型是当前值可以被自身过去p阶滞后的数据所解释，p决定了需要几个过去值来预测当前值。

$$x_{i+1}={\varphi }_1{x}_i+{\varphi }_2{x}_{i-1}+\cdots +{\varphi }_p{x}_{i-p+1}+{\xi }_{i+1}\text{\hspace{1em}(1)}$$

AR(1)模型

$$x_{i+1}=\mu +{\varphi }_1{x}_i+{\xi }_{i+1}$$

如果上面的$\varphi =1$,这个时间蓄力就是随机游走；如果$\varphi =0$，这个时间序列就是白噪声；如果$-1<\varphi <1$,这个时间序列就是平稳的。

AR(2)模型

$$x_{i+1}=\mu +{\varphi }_1{x}_i+{\varphi }_2{x}_{i-1}+{\xi }_{i+1}$$

AR(3)模型

$$x_{i+1}=\mu +{\varphi }_1{x}_i+{\varphi }_2{x}_{i-1}+{\varphi }_2{x}_{i-2}+{\xi }_{i+1}$$

模拟AR(1)模型

fig, ax = plt.subplots(ncols=1, nrows=4, figsize=(10, 10), dpi=300, sharex=True)
# AR(1) MA(1) model:AR parameter = +0.9


ar1 = np.array([1, -0.9]) # We choose -0.9 as AR parameter is +0.9
ma1 = np.array([1])
AR1 = ArmaProcess(ar1, ma1)
sim1 = AR1.generate_sample(nsample=1000)
ax[0].set_title('AR(1) model: AR parameter = +0.9',fontsize=13)
ax[0].plot(sim1)


# We will take care of MA model later
# AR(1) MA(1) AR parameter = -0.9
ar2 = np.array([1, 0.9]) # We choose +0.9 as AR parameter is -0.9
ma2 = np.array([1])
AR2 = ArmaProcess(ar2, ma2)
sim2 = AR2.generate_sample(nsample=1000)
ax[1].set_title('AR(1) model: AR parameter = -0.9',fontsize=13)
ax[1].plot(sim2)


# AR(2) MA(1) AR parameter = 0.9
plt.subplot(4,1,3)
ar3 = np.array([2, -0.9]) # We choose -0.9 as AR parameter is +0.9
ma3 = np.array([1])
AR3 = ArmaProcess(ar3, ma3)
sim3 = AR3.generate_sample(nsample=1000)
ax[2].set_title('AR(2) model: AR parameter = +0.9',fontsize=13)
ax[2].plot(sim3)



# AR(2) MA(1) AR parameter = -0.9
plt.subplot(4,1,4)
ar4 = np.array([2, 0.9]) # We choose +0.9 as AR parameter is -0.9
ma4 = np.array([1])
AR4 = ArmaProcess(ar4, ma4)
sim4 = AR4.generate_sample(nsample=1000)
ax[3].set_title('AR(2) model: AR parameter = -0.9',fontsize=13)
ax[3].plot(sim4)

对模拟数据做预测

model = ARMA(sim1, order=(1,0))
result = model.fit()
print(result.summary())
print("μ={} ,ϕ={}".format(result.params[0],result.params[1]))

可以看出来，上面的ϕ还是非常接近我们设置的模拟值的

# 将预测的数据画出来
fig, ax = plt.subplots(figsize=(10, 4), dpi=300)
# Predicting simulated AR(1) model
result.plot_predict(start=900, end=1010,ax=ax)
plt.show()

# 查看预测的rmse值
rmse = math.sqrt(mean_squared_error(sim1[900:1011], result.predict(start=900,end=999)))
print("The root mean squared error is {}.".format(rmse))

fig, ax = plt.subplots(figsize=(10, 3), dpi=300)

# 预测humidity level of Montreal
humid = ARMA(humidity["Montreal"].diff().iloc[1:].values, order=(1,0))
res = humid.fit()
res.plot_predict(start=1000, end=1100, ax=ax)
plt.show()


# 通过rmse查看预测效果
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))

上面的结果还不错，我们再来继续预测谷歌的数据

fig, ax = plt.subplots(figsize=(10, 3), dpi=300)

# Predicting closing prices of google
humid = ARMA(google["Close"].diff().iloc[1:].values, order=(1,0))
res = humid.fit()
res.plot_predict(start=900, end=1010, ax=ax)
plt.show()

MA模型

作为自回归模型的替代，q阶移动平均模型方程MA(q)中，左边的$x_t$是由右侧的白噪声$w_t$线性组合得到的。

$$x_t=w_t+{\theta }_1{w}_{t-1}+{\theta }_2{w}_{t-2}+{\theta }_3{w}_{t-3}+\cdots +{\theta }_q{w}_{t-q}$$

MA(1)模型

$$x_t=w_t+{\theta }_1{w}_{t-1}$$
上面的公式表示，今天的值就是今天的噪声加上昨天的一个噪声乘以一个$\theta$。

模拟MA(1)模型

fig, ax = plt.subplots(figsize=(10, 3), dpi=300)
ar1 = np.array([1])
ma1 = np.array([1, -0.5])
MA1 = ArmaProcess(ar1, ma1)
sim1 = MA1.generate_sample(nsample=1000)
ax.plot(sim1)
plt.show()

对模拟的MA数据做预测

model = ARMA(sim1, order=(0,1))
result = model.fit()
print(result.summary())
print("μ={} ,θ={}".format(result.params[0],result.params[1]))

# Forecasting and predicting montreal humidity
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]))
fig, ax = plt.subplots(figsize=(10, 3), dpi=300)

result.plot_predict(start=1000, end=1100, ax=ax)
plt.show()

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))

ARMA模型

我们选择处理的是平稳时间序列中的自回归、移动平均和混合自回归移动平均(ARMA)模型的一般情况。

ARMA(p,q)模型公式如下：

$$x_{i+1}={\varphi }_1{x}_i+{\varphi }_2{x}_{i-1}+\cdots +{\varphi }_p{x}_{i-p+1}+w_t+{\theta }_1{w}_{t-1}+{\theta }_2{w}_{t-2}+{\theta }_3{w}_{t-3}+\cdots +{\theta }_q{w}_{t-q}$$

使用ARMA模型做预测

# Forecasting and predicting microsoft stocks volume
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.params[2]))
fig, ax = plt.subplots(figsize=(10, 3), dpi=300)

result.plot_predict(start=1000, end=1100,ax=ax)
plt.show()


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))

可以看出来，ARMA模型的效果比AR或者MA模型都要好的多

ARIMA模型

fig, ax = plt.subplots(figsize=(10, 3), dpi=300)

# Predicting the microsoft stocks volume
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, ax=ax)
plt.show()


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))

Var模型

向量自回归模型

# Predicting closing price of Google and microsoft
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, ax = plt.subplots(figsize=(10, 10), dpi=300)
fig = plt.figure(figsize=(16,10), dpi=300)
result.plot_diagnostics(fig=fig)
# calculating error
rmse = math.sqrt(mean_squared_error(train_sample.iloc[1:1002].values, predicted_result.values))
print("The root mean squared error is {}.".format(rmse))

SARIMA模型

时间序列季节自回归求和滑动平均（SARIMA）

# Predicting closing price of Google'
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=False)
print(result.summary())
predicted_result = result.predict(start=0, end=500)
fig = plt.figure(figsize=(16,10), dpi=300)
result.plot_diagnostics(fig=fig)
# calculating error
rmse = math.sqrt(mean_squared_error(train_sample[1:502], predicted_result))
print("The root mean squared error is {}.".format(rmse))

fig, ax = plt.subplots(figsize=(10, 3), dpi=300)
ax.plot(train_sample[1:502],color='red')
ax.plot(predicted_result,color='blue')
ax.legend(['Actual','Predicted'])
fig.suptitle('Google Closing prices')
plt.show()

未观察到的成分模型(又称结构时间模型)

# Predicting closing price of Google'
train_sample = google["Close"].diff().iloc[1:].values
model = sm.tsa.UnobservedComponents(train_sample,'local level')
result = model.fit(maxiter=1000,disp=False)
print(result.summary())
predicted_result = result.predict(start=0, end=500)
fig = plt.figure(figsize=(16,10), dpi=300)

result.plot_diagnostics(fig=fig)
# calculating error
rmse = math.sqrt(mean_squared_error(train_sample[1:502], predicted_result))
print("The root mean squared error is {}.".format(rmse))

fig, ax = plt.subplots(figsize=(10, 3), dpi=300)
ax.plot(train_sample[1:502],color='red')
ax.plot(predicted_result,color='blue')
ax.legend(['Actual','Predicted'])
fig.suptitle('Google Closing prices')
plt.show()

动态因子模型

# Predicting closing price of Google and microsoft
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=False)
print(result.summary())
predicted_result = result.predict(start=0, end=1000)
fig = plt.figure(figsize=(16,10), dpi=300)
result.plot_diagnostics(fig=fig)
# calculating error
rmse = math.sqrt(mean_squared_error(train_sample.iloc[1:1002].values, predicted_result.values))
print("The root mean squared error is {}.".format(rmse))

说明

本文文字参考的是：https://www.kaggle.com/thebrownviking20/everything-you-can-do-with-a-time-series#4.-Modelling-using-statstools

代码和数据经过整理，存放在：https://gitee.com/yuanzhoulvpi/time_series

阅读更多