Python与时间序列

时间序列-统计知识

1. 什么是时间序列。
2. 在python里面如何导入时间序列数据。
3. 什么是面板数据。
4. 可视化时间序列（包括时间序列的区域填充图，季节性时间序列图，箱型图等）。
5. 时间序列的几种模式以及模式的分解。
6. 时间序列的平稳性的介绍、原理、解释等；以及让时间序列平稳的方法。
7. 白噪声和平稳时间序列的区别。
8. 时间序列的趋势性、季节性如何提取、去除。
9. 如何处理带有缺失值的时间序列数据。
10. 时间序列的自相关、偏自相关数据。
11. 时间序列的滞后图。
12. 时间序列如何平滑。

1. 什么是时间序列

同一统计指标数值按照时间先后顺序排列而成的数据。本质上是反映一个变量随时间序列变化的趋势。

1. 简单的例子就像是学生每一年的身高数据，按照时间顺序记录📝下来的数据也是一个时间序列。变量是我们的身高;
2. 我们支付宝或者微信的零钱，每一天每一个月都有一个实际的值。按照时间顺序记录下来，也是一个时间序列（没做过，每一次打开支付宝都是在流泪）。

2. 在python里面如何导入时间序列数据

一般python导入csv、excel格式的文件都是可以使用pandas包导入，如果时间格式需要被处理，可以使用datetime包来处理。

比如这样：

import numpy as np
import pandas as pd
import datetime

df = pd.read_csv('../datasets/a10.csv')
df['date'] = df['date'].apply(lambda x: datetime.datetime.strptime(x, "%Y-%m-%d"))
df.head()

输出的结果如下：

3. 什么是面板数据

有多个时间节点，且每一个时间节点上有多个维度的数据。

举一个简单的例子，如果你是一个学生，每天都要考试，2021-01-01日你语文90、数学100、英语102；2021-01-02日你语文10、数学110、英语122 ……这样的就是就是一个面板数据，有多个时间节点，且每一个时间节点，有多个维度（语文、数学、英语）；

举一个现实的例子；后来这个学生长大了，每天有自己的经济来源了，他还是喜欢记录数据，他会每天都会这么记录：2021-01-01日给女友花费1000元、兼职20元、日常开销200元；2021-01-02日给女友10元、兼职200元、日常开销10元；…… 这也是一个面板数据。（怎么有点心酸呢）。

4. 可视化时间序列

1. 最常见的就是时序图：

import matplotlib.pyplot as plt

df = pd.read_csv('../datasets/a10.csv')
df['date'] = df['date'].apply(lambda x: datetime.datetime.strptime(x, "%Y-%m-%d"))

def plot_df(df, x, y, title="", xlabel="Date", ylabel="Value", dpi=100):
    plt.figure(figsize=(16, 5), dpi=dpi)
    plt.plot(df[x], df[y], color='tab:red')
    plt.gca().set(title=title, xlabel=xlabel, ylabel=ylabel)
    plt.show()
plot_df(df, x="date", y="value")

2. 区域填充图：

fig, ax = plt.subplots(1, 1, figsize=(12, 5), dpi=100)
ax.fill_between(x=df['date'], y1=df['value'], y2=-df['value'], alpha=0.5, linewidth=2, color='seagreen')
ax.set_title("demo data", fontsize=16)
ax.hlines(y=0, xmin=np.min(df['date']), xmax=np.max(df['date']), linewidth=0.5)
plt.show()

3. 季节性图

df['year'] = df['date'].apply(lambda x: x.year)
df['month'] = df['date'].apply(lambda x: int(x.strftime('%m')))

fig, ax = plt.subplots(figsize=(10, 6))
for index, tempyear in enumerate(df['year'].unique().tolist()):
    tempdata = df.loc[df['year'] == tempyear].sort_values(by=['date'])
    ax.plot(tempdata['month'], tempdata['value'], label=tempyear)
    ax.text(tempdata.tail(1)['month'], tempdata.tail(1)['value'], tempyear)

font = {'family': 'serif',
        'color': 'darkred',
        'weight': 'normal',
        'size': 16,
        }
ax.set_ylabel("$Drug Sales$", fontdict=font)
ax.set_xlabel('$Month$', fontdict=font)
ax.set_title("Seasibak Plot of Drug Sales Time Series",
             fontdict=font)
x_axis_ticks = np.arange(start=1, stop=13)
ax.set_xticks(x_axis_ticks)
ax.set_xticklabels(x_axis_ticks)

4. 对季节性做下钻，可以使用boxplot

fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(16, 5))
sns.boxplot(x='year', y='value', data=df, ax=ax[0])
sns.boxplot(x='month', y='value', data=df, ax=ax[1])

ax[0].set_title("Year-wise BoxPlot \n(The Trend)", fontsize=18)
xlabel = ax[0].get_xticklabels()
ax[0].set_xticklabels(xlabel, rotation=45)
ax[1].set_title("Month-wise BoxPlot \n(The Trend)", fontsize=18)

5. 时间序列中的模式

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

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

但是不是说所有的时间序列都必须具备这3个部分。时间序列可能没有明显的趋势、可能没有明显的周期性。或者两个都没有。因此，可以将时间序列看成趋势、周期性、误差项的组合。

1. 几种模式的对比

fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(20, 4), dpi=100)
pd.read_csv("../datasets/guinearice.csv", parse_dates=['date'], index_col=['date']).plot(title="Trend Only", legend=False,ax=ax[0])
pd.read_csv("../datasets/sunspotarea.csv", parse_dates=['date'], index_col=['date']).plot(title="Seasonality Only",legend=False, ax=ax[1])
pd.read_csv("../datasets/AirPassengers.csv", parse_dates=['date'], index_col=['date']).plot(title="Trend and Seasonality",legend=Fal

5.1 时间序列中的模式分解

基于原始的时间序列的趋势项和季节项，时间序列模型可以被分为加法模型和乘法模型.

1. 乘法型：时间序列值 = 趋势项 * 季节项 * 误差项。
2. 加法型：时间序列值 = 趋势项 + 季节项 + 误差项。

对时间序列的分解，可以从趋势项、季节项、 误差项的乘法和加法的角度进行。代码方面，我们使用statsmodels包的seasonal_decompose函数。

from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse

# 导入数据
df = pd.read_csv("../datasets/a10.csv", parse_dates=['date'], index_col=['date'])

# 乘法模型
result_mul = seasonal_decompose(df['value'], model="multilicative", extrapolate_trend='freq')
result_add = seasonal_decompose(df['value'], model="additive", extrapolate_trend='freq')

# 画图
fig, ax = plt.subplots(ncols=2, nrows=4, figsize=(14, 7), sharex=True)


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

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

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

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


plot_decompose(result=result_add, ax=ax, index=0, title="Additive Decompose", fontdict=font)
plot_decompose(result=result_mul, ax=ax, index=1, title="Multiplicative Decompose", fontdict=font)

1. 在上面的代码中，设置extrapolate_trend='freq'是为了填充 季节项、误差项中开头的缺失值。
2. 在上图中，可以发现加法模型里面的误差项还有一部分时间序列的模式没有被提取完。但是乘法模型的误差项(或者叫残差项)基本上看不出来任何信息了，说明乘法模型对这个数据来说，有更强提取数据信息的能力。
3. 乘法模型的季节项、趋势项、误差项数据都保存在result_mul里面。我们将这几个数据提取出来。

df_reconstructed = pd.concat([result_mul.seasonal, result_mul.trend, result_mul.resid, result_mul.observed], axis=1)
df_reconstructed.columns = ['Seasonal', 'Trend', 'Resid', 'Actual_value']
df_reconstructed.head()

可以检查一下上面的数据;基本上可以确定df_reconstructed['Seasonal'] * df_reconstructed['Trend'] * df_reconstructed['Resid'] = df_reconstructed['Actual_value'], 或者我们用均方误差算一下：

value = np.sqrt(np.sum((df_reconstructed['Seasonal'] * df_reconstructed['Trend'] * df_reconstructed['Resid'] -df_reconstructed['Actual_value']) ** 2))

value

6. 时间序列的平稳性

1. 平稳性是时间序列的一个属性，一个平稳的时间序列指的是这个时间序列和时间无关，也就是说这个时间序列不是时间的一个函数。
2. 也就是说一个时间序列是平稳的，就是这个时间序列的几个统计量：均值、方差、自相关系数都是一个常数，和时间无关。

6.1 如何让一个时间序列平稳

1. 对一个时间序列做一次或者多次差分。
2. 季节性延迟（或者叫多阶差分）。
3. 取时间序列的第N个根。
4. 上面三种方法混合使用。

最方便最常用的平稳一个时间序列，就是多次使用一阶差分，直到时间序列平稳。

6.2在做预测之前为什么需要对非平稳时间序列做平稳

1. 预测一个平稳的时间序列相对容易，且结果更加可靠。
2. 一个更加重要的原因是：自回归预测模型本质上是一个线性回归模型，它利用时间序列的滞后项做为预测变量。我们都知道，在线性回归里面，如果预测变量的X都是互不相关的，那么线性回归预测的效果最好。因此对序列进行平稳化就解决了变量之间的相关性问题， 从而消除来时间序列的自相关，使得预测模型中的预测变量几乎相互独立。

现在已经知道的时间序列的平稳性的重要性，如何判断一个时间序列是否平稳呢？

6.3检验时间序列是否平稳

检验时间序列是否为平稳性时间序列可以有这几种方法：

1. 查看时间序列的时序图；
2. 将时间序列分为2个或者多个连续片段，然后计算对应的统计量（比如均值、方差、自相关系数等），如果片段的统计量都是明显不相等的话，说明肯定不是平稳性。
3. 不管怎么样，需要使用一个量化的指标去衡量一个时间序列是否平稳，可以使用单位根检验方法来对时间序列做平稳性检验，检验时间序列是否平稳并且具有单位根。

还有很多单位根检验方法的增强版本：

1. Augmented Dickey Fuller test(ADH Test)
2. Kwiatkowski-Phillips-Schmidt-Shin – KPSS test (trend stationary)
3. Philips Perron test (PP Test)

最广泛使用的方法是ADF test；原假设是时间序列有单位根并且是非平稳的。如果ADF test的P-value小于显著水平（0.05），就可以拒绝原假设。KPSS test和ADH test恰恰相反，是为了证明时间序列有一个确定的趋势。

from statsmodels.tsa.stattools import adfuller, kpss

df = pd.read_csv("../datasets/a10.csv", parse_dates=['date'])

# ADF test
result = adfuller(df.value.values, autolag='AIC')

print('*' * 50)
print(f"ADF Statistic: {result[0]}; p-value: {result[1]}")

for key, value in result[4].items():
    print('Crital Values:')
    print(f"{key}, {value}")

# KPSS Test

result = kpss(df.value.values, regression='c', nlags='auto')

print('*' * 50)
print(f'\nKPSS Stattistic: {result[0]:.10f}; p-value: {result[1]:.3f}')

for key, value in result[3].items():
    print('Crital Values:')
    print(f"{key}, {value}")

结果：

6.4白噪声和平稳时间序列的区别

和平稳性时间序列一样，白噪声时间序列的每一个值和时间也是无关的，也就是说白噪声的均值和方差不随着时间改变。唯一一个不同点就是白噪声是完全随机的、且均值为0的时间序列。

白噪声中没有任何时间序列模式，生活中遇到的白噪声就像是FM收音机中空白声音。数学上的白噪声就是一串均值为0的完全随机的随机数。

fig, ax = plt.subplots(figsize=(10, 4))
randval = np.random.randn(1000)
pd.Series(randval).plot(title="Random White Noise", color='k', ax=ax)

8.时间序列的去趋势

去趋势时间序列指的是从时间序列中去趋势成分。去除时间序列趋势有多种方法。

1. 从时间序列中获得最佳拟合线。可以使用线性回归来拟合时间序列，有时候可能还用得到二次项(x^2);
2. 将刚才我们提取的时间序列成分中趋势项部分移除；
3. 使用像是Baxter-King过滤器或者Hodrick-Prescott过滤器去除移动平均趋势线或者周期性成分。

8.1去时间序列的趋势项趋势

# use scipy: Subtract the Line of best fit
from scipy import signal

df = pd.read_csv("../datasets/a10.csv", parse_dates=['date'])
detrend = signal.detrend(df.value.values)

fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(15, 4))
ax[0].plot(detrend)
ax[0].set_title("Drug Sales detrended by subtracting the least squares fit")

# Using  statmodels: Subtracting the Trend Component
from statsmodels.tsa.seasonal import seasonal_decompose

df = pd.read_csv('../datasets/a10.csv', parse_dates=['date'], index_col='date')
result_mul = seasonal_decompose(df['value'], model="multilicative", extrapolate_trend='freq')
detrended = df.value.values / result_mul.trend
ax[1].plot(detrended)
ax[1].set_title("Drug Sales detrended by subtracting the trend component")

8.2去时间序列的季节性趋势

fig, ax = plt.subplots(figsize=(8, 4))

# Subtracing the Trend Component
from statsmodels.tsa.seasonal import seasonal_decompose

df = pd.read_csv('../datasets/a10.csv', parse_dates=['date'], index_col='date')
result_mul = seasonal_decompose(df['value'], model="multilicative", extrapolate_trend='freq')
detrended = df.value.values / result_mul.seasonal
ax.plot(detrended)
ax.set_title("Drug Sales Deseasonlized")

9.检验时间序列的季节性

最常用的方式是将时间序列的图画出来，然后从图中查看重复的模式和固定时间间隔，一般周期性都是由日期或者时间决定：

1. 一天中的小时；
2. 一个月中的每天；
3. 每周的；
4. 每月的；
5. 每年的。

如果想对季节性做更加准确的季节性检验，就使用自相关函数图(ACF plot)；但是当存在相当明显的季节性模式的时候，ACF图会在季节性窗口长度的倍数处出现明显的峰值。

比如，前面的毒品销售是每月的系列，每年都有重复的模式，比如在下面的第12，第24，第36位置，都可以看到线有个小峰值。但是必须注意：在实际数据集中，这些强累的模式几乎不会被注意到，而且可能还会被噪音干扰。

from pandas.plotting import autocorrelation_plot

df = pd.read_csv("../datasets/a10.csv")
fig, ax = plt.subplots(figsize=(10, 4))
autocorrelation_plot(df.value, ax=ax)

ybound = ax.get_ybound()
for i in 12 * np.arange(start=1, stop=13):
    ax.vlines(x=i, ymin=ybound[0], ymax=ybound[1], colors='red', linestyles='--', linewidth=0.6)
    ax.text(x=i, y=ybound[1], s=f"{i:d}th", fontdict={'size': 8})

ax.grid(False)

相应的，如果需要统计检验，可以使用CH-test来决定是否需要对时间序列做季节差分。

10.对带有缺失值的时间序列处理

现实中的时间序列数据可能会缺少时间（比如少了几天，少了几年），这也就意味着，这些时间段内的数据没有被捕捉到或者不能用，在这些缺失的时间段范围内，测量值可能为0，这样的情况下，可以使0来填充这个时间段内的数据。对时间序列做填充的时候，不能简单的对缺失值用均值填充，应该使用向前填充，用后一个数据去填前面的数据。处理现实的数据，需要做很多尝试，来确定哪一个填充效果最好，比如有：

1. 向前\向后填充；
2. 线性插值、二次插值、最近邻插值、使用季节性插值。

为了衡量插值的效果，做了一次小实验，手动对数据做缺失值处理，然后使用不同方法做插值。

# 生成随机数
from scipy.interpolate import interp1d
from sklearn.metrics import mean_squared_error

df_orig = pd.read_csv("../datasets/a10.csv", parse_dates=['date'], index_col=['date'])

# 制作缺失值数据
np.random.seed(2021)
df = df_orig.copy()
df.iloc[np.random.randint(low=0, high=df_orig.shape[0], size=30)] = None

# 画图
fig, ax = plt.subplots(nrows=7, ncols=1, figsize=(10, 12), sharex=True)

# 1. 原始数据 -------------------------------
df_orig.plot(title="Acutal", ax=ax[0], label='Acutal', color='red', style='.-')
df.plot(ax=ax[0], label='Acutal', color='green', style='.-')
ax[0].legend(['Missing Data', 'Available Data'])

# 2. 向前填充 -------------------------------
df_ffill = df.ffill()
error = mean_squared_error(df_orig['value'], df_ffill['value'])
df_ffill['value'].plot(title=f"Forward Fill (MSE: {error:.3f})", ax=ax[1])

# 3. 向后填充 -------------------------------
df_bfill = df.bfill()
error = mean_squared_error(df_orig['value'], df_bfill['value'])
df_ffill['value'].plot(title=f"Forward Fill (MSE: {error:.3f})", ax=ax[2])

# 4. 线性插值 -------------------------------
df['rownum'] = np.arange(df.shape[0])
df_nana = df.dropna(subset=['value'])

f = interp1d(df_nana['rownum'], df_nana['value'])
df['linear_fill'] = f(df['rownum'])

error = mean_squared_error(df_orig['value'], df['linear_fill'])
df['linear_fill'].plot(title=f"Linear Fill (MSE: {error:.3f})", ax=ax[3])

# 5. 三次方插值 -------------------------------
# df['rownum'] = np.arange(df.shape[0])
# df_nana = df.dropna(subset=['value'])
f = interp1d(df_nana['rownum'], df_nana['value'], kind='cubic')
df['cubic_fill'] = f(df['rownum'])

error = mean_squared_error(df_orig['value'], df['cubic_fill'])
df['cubic_fill'].plot(title=f"cubic Fill (MSE: {error:.3f})", ax=ax[4])


# 6. 邻近均值插值 -------------------------------
def knn_mean(ts, n):
    out = np.copy(ts)

    for i, val in enumerate(ts):
        if np.isnan(val):
            n_by_2 = np.ceil(n / 2)
            lower = np.max([0, int(i - n_by_2)])
            upper = np.min([len(ts) + 1, int(i + n_by_2)])
            ts_near = np.concatenate([ts[lower:i], ts[i:upper]])
            out[i] = np.nanmean(ts_near)

    return out


df['knn_mean'] = knn_mean(df.value.values, 8)
error = mean_squared_error(df_orig['value'], df['knn_mean'])
df['knn_mean'].plot(title=f"KNN Mean (MSE: {error:.3f})", ax=ax[5])


# 7.季节性均值插值 -------------------------------

def seasonal_mean(ts, n, lr=0.7):
    """
    计算时间序列每一个季节性的均值，填充到季节性里面的缺失值。
    ts：一个时间序列；
    n: 季节性的周期
    """
    out = np.copy(ts)
    for i, val in enumerate(ts):
        if np.isnan(val):
            ts_seas = ts[i - 1::-n]
            if np.isnan(np.nanmean(ts_seas)):
                ts_seas = np.concatenate(ts[i - 1::-n], ts[i::n])
            out[i] = np.nanmean(ts_seas) * lr

    return out


df['seasonal_mean'] = seasonal_mean(df.value, n=12, lr=1.25)
error = mean_squared_error(df_orig['value'], df['seasonal_mean'])
df['seasonal_mean'].plot(title=f"Seasonal Mean (MSE: {error:.3f})", ax=ax[6])

填充缺失值还有很多方法，填充的效果是取决于你对填充精度的要求。

1. 如果还有很多辅助变量，可以使用随机森林或者k邻近等模型做预测填充；
2. 可以使用过去的数据填充未来的；可以使用现在的数据去填充过去的。
3. 利用时间序列的季节性填充数据；
4. 利用时间序列的周期性对数据u哦填充。

11.什么是自相关和偏自相关

1. 自相关是时间序列和它自己滞后n阶时序数据的相关性。如果一个时间序列有显著的自相关性，说明使用历史数据做预测效果会非常好。
2. 偏自相关性是和自相关系数差不多，但是不包括中间滞后的相关贡献。

from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

df = pd.read_csv('../datasets/a10.csv')

# Calculate ACF and PACF upto 50 lags
# acf_50 = acf(df.value, nlags=50)
# pacf_50 = pacf(df.value, nlags=50)

# Draw Plot
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(16, 3), dpi=100)
plot_acf(df.value.tolist(), lags=50, ax=ax[0])
ax[0].set_title("acf")
plot_pacf(df.value.tolist(), lags=50, ax=ax[1])
ax[1].set_title("pacf")

11.1偏自相关系数到底是什么

有一个时间序列叫,他的一阶滞后叫，他的二阶滞后叫，他的三阶滞后叫。

这里有一个公式：.

那么这里就是第三阶偏自相关系数。

12.滞后项画图

滞后项画图，就是一个时间序列和他自己滞后n阶的散点图，通常可以用来检查自相关性，如果可以从图中看到一些相关性，说明数据有一点相关性的。如果从滞后项图看不出来什么数据信息，很有可能说明数据是白噪声。

from pandas.plotting import lag_plot

# import data
ss = pd.read_csv("../datasets/sunspotarea.csv")
a10 = pd.read_csv("../datasets/a10.csv")

# plot Sun Spot Area
fig, ax = plt.subplots(ncols=4, nrows=1, figsize=(10, 3), sharex=True, sharey=True, dpi=100)

for i, ax in enumerate(ax.flatten()):
    lag_plot(ss['value'], lag=i + 1, ax=ax, c='firebrick')
    ax.set_title(f"Lag ({i + 1})")
fig.suptitle("Lag Plot of Sun Spot Area\n(Points get wide and scattered with increasing lag -> lesser correlation)",
             fontdict=font, y=1.12)

# plot  Drugs Sales
fig, ax = plt.subplots(ncols=4, nrows=1, figsize=(10, 3), sharex=True, sharey=True, dpi=100)

for i, ax in enumerate(ax.flatten()):
    lag_plot(a10['value'], lag=i + 1, ax=ax, c='firebrick')
    ax.set_title(f"Lag ({i + 1})")
fig.suptitle("Lag Plots of Drug Sales", fontdict=font, y=1.05)

从下图可以观察到，sun spot的滞后阶数远小于季节性，所以滞后1、2、3、4的数据和原始数据看起来没什么关系，可以说明时间序列前后可能没有什么相关性。但是下面的毒品售卖滞后1、2、3、4阶和原始数据看起来有很强的相关性，需要向下继续研究。

13.如何对时间序列做平滑

对时间序列做平滑的优点：

1. 减少噪声对时间序列的影响；获得除去噪声后的时间序列数据；
2. 时间序列平滑后的数据可以用来解释原序列的一些特征；
3. 更好的可视化潜在的趋势。

如何平滑时间序列：

1. 移动平均；
2. 做一个局部平滑（LOESS）；
3. 做一个局部加权回归（LOWESS ）。

做移动平均的时候，需要注意窗口的大小。如果窗口过大，导致时间序列过于平滑，以至于还消除了一些其他信息。

LOESS是局部回归的意思，取局部数据做回归；LOWESS是局部加权回归。

from statsmodels.nonparametric.smoothers_lowess import lowess

# 导入数据
df_orig = pd.read_csv("../datasets/elecequip.csv", parse_dates=['date'], index_col=['date'])

# 1. 移动平均
df_ma = df_orig.value.rolling(3, center=True, closed='both').mean()

# 2. loess 平滑 (5% and 15%)
# df_loess_5 = pd.DataFrame(lowess(df_orig.value, np.arange(len(df_orig.value))[:,1]), index=df_orig.index, columns=['value'])
df_loess_5 = pd.DataFrame(lowess(endog=df_orig.value, exog=np.arange(len(df_orig.value)), frac=0.05)[:,1],index=df_orig.index, columns=['value'])
df_loess_15 = pd.DataFrame(lowess(endog=df_orig.value, exog=np.arange(len(df_orig.value)), frac=0.15)[:,1],index=df_orig.index, columns=['value'])

# plot
fig, ax = plt.subplots(4,1, figsize=(7, 7), sharex=True, dpi=100)
df_orig['value'].plot(ax=ax[0], color='k', title='Original Series')
df_loess_5['value'].plot(ax=ax[1], title='Loess Smoothed 5%')
df_loess_15['value'].plot(ax=ax[2], title='Loess Smoothed 15%')
df_ma.plot(ax=ax[3], title='Moving Average (3)')
fig.suptitle('How to Smoothen a Time Series', y=0.95, fontsize=14)

时间序列-相关系数

1. 时间序列中最基本的概念；
2. 什么叫lag（滞后）；
3. 时间序列的公式；
4. 统计基础知识；
5. 自相关系数（公式、复现、调包、可视化）；
6. 偏自相关系数（公式、复现、调包、可视化）；

1.时间序列中最基本的概念

很多人看不懂时间序列的公式、很多人看不懂时间序列的数据、代码、数学函数。本文就是将数据、公式、代码、数学函数串联起来，让你不仅仅了解抽象的数学公式，也能将抽象的内容联系到具体的数据和python代码中，并且也不只是单纯的调用python包，而是自己手写算法。

下面的代码是近20年城镇就业人员数据(万人)，数据来自中国统计局。

import numpy as np
import pandas as pd
import datetime
import matplotlib.pyplot as plt

data = pd.read_csv("../datasets/近20年城镇就业人员数据(万人).csv", index_col=['year'])
data.head(4)

# 可视化
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(data.index, data['value'], 'o-')
ax.set_title(u"Data on urban employees in the past 20 years")
ax.set_ylabel("Ten thousand people")
ax.set_xlabel("year")
ax.set_xticks(data.index)
ax.set_xticklabels(data.index, rotation=45)

2.什么叫lag（滞后）

经常在时间序列里面可以看到lag这个东西，大白话就是叫滞后；举个例子来说：

1. lag1代表滞后数据1项，也就是2002年的数据滞后1年得到的就是2001年的数据；2003年的数据滞后1年就是2002年的数据；
2. lag2代表滞后数据2项，也就是2003年的数据滞后2年得到2001年的数据；2004年的数据滞后2年就是2002年的数据；
3. 依次类推，滞后数据就是这么计算的。

下面使用python对上面的【近20年城镇就业人员数据(万人)】数据做个滞后处理。

# 使用pandas的shift可以做lag
data1 = data.copy()
# data1['value_lag1'] = data1['value'].shift(1)
# data1['value_lag2'] =  data1['value'].shift(2)
for i in [1, 2, 3, 4, 5]:
    data1[f"value_lag{i}"] = data1['value'].shift(i)
data1.head(9)

数据解读：

1. 因为2001年之前是没有数据了，所以2001年滞后一项得到的数据是空数据；
2. 因为2002年之前是没有数据，所以2002年滞后两项也是没有数据的。
3. 依次类推。

3.时间序列的公式

经常看到时间序列的公式写的乱七八糟的，这是一堆公式，那是一堆公式，活活把人绕晕。现在为了我们说明需要，我们定义几个简单的公式：

假设我们的时间序列是value = [4,3,2,10,3]。我们把这个时间序列定义为 Xt；那么这个时间序列里面每一个元素对应的符号依次是：X1,X2,X3等。

4.统计基础知识

5.自相关系数

5.1自相关基础计算

在了解上面的【统计基础知识】后，再理解自相关性基本上就不难了，自相关是变量与自身在不同时间滞后下的相关性。

对于延迟k阶的时间序列。自协方差为1，自相关性系数为2

5.2调用包计算acf数值

# 使用statsmodels包的acf函数
from statsmodels.tsa.stattools import acf

5.3自己根据公式写函数

def cal_acf(x, nlags):
    """
    按照公式自己写个acf函数
    只是用来帮助理解，不建议用在别的环境中
    """
    x = np.array(x)
    mean_x = np.mean(x)
    length_x = x.shape[0]
    c_0 = np.mean((x-mean_x) **2)
    c_k = np.sum((x[:(length_x-nlags)] - mean_x) * (x[nlags:] - mean_x)) / length_x
    r_k = c_k / c_0
    return r_k

5.4两种结果对比

# 结果汇总
pd.DataFrame({'index':np.arange(11),
             'value_by_myself':[cal_acf(x=data['value'], nlags=i) for i in range(11)],
             'value_by_statsmodels':acf(data.value,nlags=10)})

5.5自相关系数画图

画图有两种方法：

1. 直接使用statsmodels的plot_acf;
2. 使用acf函数来计算得到数值，得到数值后再进行可视化，这种一般便于自己定制，也方便自己做更多的细节调整之类的。

# 自己画一个
# 自己计算
acf_value, acf_interval, _, _ = acf(data.value,nlags=14,qstat=True,alpha=0.05, fft=False)

xlabel = np.arange(start=0, stop=acf_value.shape[0], dtype='float')

fig, ax = plt.subplots(nrows=2, figsize=(10,6), sharex=True, dpi=100)
ax[0].hlines(y=0, xmin=np.min(xlabel)-2, xmax=np.max(xlabel)+2, colors='red', linewidth=0.5)
ax[0].scatter(x=xlabel, y=acf_value, c='red')
ax[0].vlines(x = xlabel, ymin=0, ymax=acf_value, colors='red', linewidth=0.5)
xlabel[1] -= 0.5
xlabel[-1] += 0.5
ax[0].fill_between(x=xlabel[1:], y1=acf_interval[1:,0] - acf_value[1:], y2=acf_interval[1:, 1]-acf_value[1:], alpha=0.25, linewidth=0, color='red')
ax[0].set_title("acf plot by myself")


# 使用别人写好的函数

from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data['value'], ax=ax[1])
ax[1].set_title("acf plot by statsmodels")
ax[1].set_xlim(-1, np.max(xlabel)+1)

6.偏自相关系数

偏自相关系数可以被看成一种条件相关，两个变量之间有相关性是基于一系列的其他的条件变量。想象一个这样的问题：

6.1偏自相关计算公式：

6.2偏自相关计算：

1. 使用statsmodels的pacf函数，这个函数可以选择不同的计算方法，有使用ols方法的，也有使用Yule-Walker方法的。
2. 除了方法1调用包，我们还可以使用公式方法，分别使用ols方法、使用Yule-Walker方法来计算pacf的系数。实际上算出来的结果和statsmodels的pacf函数算出来的结果是一模一样的。

在下面的代码块中，我们使用statsmodels包里面的pacf、pacf_ols、pacf_yw函数计算偏自相关系数，另外我们自己通过ols理论和Yule-Walker理论计算偏自相关系数。

# 使用statsmodels包的pacf函数
from statsmodels.tsa.stattools import pacf,pacf_ols,pacf_yw


# 自己实现用ols求解pacf的函数
from numpy.dual import lstsq
from statsmodels.tools import add_constant
from statsmodels.tsa.tsatools import lagmat

def cal_my_pacf_ols(x, nlags=5):
    """
    自己实现pacf，原理使用的就是ols(最小二乘法)
    :param x:
    :param nlags:
    :return:
    """
    pacf = np.empty(nlags+1) * 0
    pacf[0] = 1.0

    xlags, x0 = lagmat(x, nlags, original="sep")
    xlags = add_constant(xlags)

    for k in range(1, nlags+1):
        params = lstsq(xlags[k:, :(k+1)], x0[k:], rcond=None)[0]
        pacf[k] = params[-1]

    return pacf




def cal_my_yule_walker(x, nlags=5):
    """
    自己实现yule_walker理论
    :param x:
    :param nlags:
    :return:
    """
    x = np.array(x, dtype=np.float64)
    x -= x.mean()
    n = x.shape[0]

    r = np.zeros(shape=nlags+1, dtype=np.float64)
    r[0] = (x ** 2).sum()/n

    for k in range(1, nlags+1):
        r[k] = (x[0:-k] * x[k:]).sum() / (n-k*1)

    from scipy.linalg import toeplitz
    R = toeplitz(c=r[:-1])
    result = np.linalg.solve(R, r[1:])
    return result

def cal_my_pacf_yw(x, nlags=5):
    """
    自己通过yule_walker方法求出pacf的值
    :param x:
    :param nlags:
    :return:
    """
    pacf = np.empty(nlags+1) * 0
    pacf[0] = 1.0
    for k in range(1, nlags+1):
        pacf[k] = cal_my_yule_walker(x,nlags=k)[-1]

    return pacf

cal_my_pacf_yw(x=data.values)

# 结果比较

pd.DataFrame({'pacf':pacf(data.value,nlags=5, method='ols'),
              'pacf_ols':pacf_ols(data.value,nlags=5),
              'pacf_yw':pacf_yw(data.value,nlags=5),
              'pacf_ols_bymyself':cal_my_pacf_ols(x=data.values),
              'pacf_yw_bymyself':cal_my_pacf_yw(x=data.values)})

效果对比：

6.3偏自相关系数画图

画图有两种方法：

1. 直接使用statsmodels的plot_pacf;
2. 使用pacf函数来计算数据，得到数据后再进行可视化，这种一般便于自己定制，也方便自己做更多的细节调整之类的。

#%%
sunspotarea = pd.read_csv("../datasets/sunspotarea.csv",parse_dates=['date'])


# 自己画一个
# 自己计算
pacf_value, pacf_interval = pacf(sunspotarea.value,nlags=15,alpha=0.05)

xlabel = np.arange(start=0, stop=pacf_value.shape[0], dtype='float')

fig, ax = plt.subplots(nrows=2, figsize=(10,6), sharex=True, dpi=100)
ax[0].hlines(y=0, xmin=np.min(xlabel)-2, xmax=np.max(xlabel)+2, colors='red', linewidth=0.5)
ax[0].scatter(x=xlabel, y=pacf_value, c='red')
ax[0].vlines(x = xlabel, ymin=0, ymax=pacf_value, colors='red', linewidth=0.5)
xlabel[1] -= 0.5
xlabel[-1] += 0.5
ax[0].fill_between(x=xlabel[1:], y1=pacf_interval[1:,0] - pacf_value[1:], y2=pacf_interval[1:, 1]-pacf_value[1:], alpha=0.25, linewidth=0, color='red')
ax[0].set_title("pacf plot by myself")


# 使用别人写好的函数


from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(sunspotarea.value, ax=ax[1], lags=15)
ax[1].set_title("pacf plot by statsmodels")
ax[1].set_xlim(-1, np.max(xlabel)+1)

Yule-Walker方程

使用Yule-Walker方法计算AR(p)

1.转置

1.1 p=1

1.2 p=2

2. Yule-Waler公式

这是一个AR(p)公式：

2.1 lag=1（滞后1）

2.2 lag=2（滞后2）

2.3 lag=k(滞后k)

2.4 lag=p(滞后p)

2.5 将上面的公式放一起

3. Yule-Walker公式求解过程

4. python计算Yule-Walker公式

代码如下：

from scipy.linalg import toeplitz
import numpy as np


def cal_my_yule_walker(x, nlags=5):
    """
    自己实现yule_walker理论
    :param x:
    :param nlags:
    :return:
    """
    x = np.array(x, dtype=np.float64)
    x -= x.mean()
    n = x.shape[0]

    r = np.zeros(shape=nlags+1, dtype=np.float64)
    r[0] = (x ** 2).sum()/n

    for k in range(1, nlags+1):
        r[k] = (x[0:-k] * x[k:]).sum() / (n-k*1)


    R = toeplitz(c=r[:-1])
    result = np.linalg.solve(R, r[1:])
    return result

def cal_my_pacf_yw(x, nlags=5):
    """
    自己通过yule_walker方法求出pacf的值
    :param x:
    :param nlags:
    :return:
    """
    pacf = np.empty(nlags+1) * 0
    pacf[0] = 1.0
    for k in range(1, nlags+1):
        pacf[k] = cal_my_yule_walker(x,nlags=k)[-1]

    return pacf


# 测试函数
data_x = np.random.randint(low=10, high=20, size=20)

# 使用yelu_walker方法计算pacf
cal_my_pacf_yw(data_x)

时间序列模型

• 介绍时间序列

• • 导入时间序列数据
  • 处理(清洗)数据
  • 对数据进行可视化
  • 时间戳和时间段
  • 使用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(4)

#%%

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

处理(清洗)数据

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

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

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

对数据进行可视化

# 将时间序列转换成特定的频率，下面的.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决定了需要几个过去值来预测当前值。

模拟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模型

模拟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模型做预测

# 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://gitee.com/yuanzhoulvpi/time_series