[人工智能]定量预测方法总结及案例实践

1 前序

常见的两类方法：

• 基于相关原则的回归分析法
• 基于惯性原则的时间序列法

回归分析的基本步骤：

• （1）重点考察一个特定的变量（因变量），而把其他变量（自变量）看作是影响这一变量的因素，并通过适当的数学模型将变量间的关系表达出来；
• （2）利用样本数据建立模型的估计方程；
• （3）对模型进行显著性检验；
• （4）通过一个或几个自变量的取值来估计或预测因变量的取值。

预测方法的选择（仅供参考）：

预测方法的评估：一种预测方法的好坏取决于预测误差的大小。预测误差是预测值与实际值的差距，度量方法有：

• 平均误差（Mean Error）
• 平均绝对误差（Mean Absolute Deviation）
• 均方误差（Mean Square Error,MSE）(常用)
• 平均百分比误差（Mean Percentage Error）
• 平均绝对百分比误差（Mean Absolute Percentage Error）

2 预测方法及案例

2.1 回归分析

2.1.1 含有哑变量的线性回归分析案例

为研究员工月工资收入与工作年限和性别之间的关系，从某公司职员中随机抽取男女各4名，他们的月工资收入与工作年限和性别之间的关系表如下：

令$y$表示月工资收入，$x_1$表示工作年限，$x_2$表示性别，性别作为哑变量引入时，回归方程如下：$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2$，于是我们可以得到：

• 女（$x_2=0$）：$y_{女性}={\beta }_0+{\beta }_1{x}_1$
• 男($x_2=1$)：$y_{男性}=({\beta }_0+{\beta }_2)+{\beta }_1{x}_1$

其中各参数的含义如下：

• ${\beta }_0$的含义是女性职工的基本月工资收入
• $({\beta }_0+{\beta }_2)$的含义是男性职工的基本月工资收入
• ${\beta }_1$的含义是工作年限每增加1年，男性或女性工资的平均增加值
• ${\beta }_2$的含义是男性职工的月工资收入与女性职工的月工资收入之间的差值，即$y_{男性}?y_{女性}=({\beta }_0+{\beta }_2)+{\beta }_1{x}_1?{\beta }_0+{\beta }_1{x}_1={\beta }_2$

python实现代码如下：

import pandas as pd
import numpy as np
import statsmodels.api as sm

data = pd.DataFrame({
    '月工资收入':[2900,3000,4800,1800,2900,4900,4200,4800],
    '工作年限':[2,6,8,3,2,7,9,8],
    '性别':['男','女','男','女','男','男','女','女']
})

# 哑变量处理
dummy_variables = pd.get_dummies(data=data['性别'].values)

X = np.column_stack(tup=(data['工作年限'].values,dummy_variables.values))
X = sm.add_constant(data=X)     # 加一列常数项
y = data['月工资收入'].values

# 用最小二乘法拟合回归方程
linear_model = sm.OLS(endog=y,exog=X)
ols_result = linear_model.fit()

# 输出拟合结果
print(ols_result.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.901
Model:                            OLS   Adj. R-squared:                  0.862
Method:                 Least Squares   F-statistic:                     22.78
Date:                Sun, 22 May 2022   Prob (F-statistic):            0.00307
Time:                        17:02:55   Log-Likelihood:                -58.036
No. Observations:                   8   AIC:                             122.1
Df Residuals:                       5   BIC:                             122.3
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        950.7246    247.685      3.838      0.012     314.029    1587.420
x1           397.5845     60.183      6.606      0.001     242.879     552.290
x2           -85.0242    231.138     -0.368      0.728    -679.184     509.135
x3          1035.7488    172.207      6.015      0.002     593.076    1478.421
==============================================================================
Omnibus:                        4.593   Durbin-Watson:                   1.536
Prob(Omnibus):                  0.101   Jarque-Bera (JB):                1.483
Skew:                           1.049   Prob(JB):                        0.477
Kurtosis:                       3.219   Cond. No.                     7.55e+16
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The smallest eigenvalue is 5.64e-32. This might indicate that there are
strong multicollinearity problems or that the design matrix is singular.

从上面结果可以看到，拟合优度R方为0.901，调整R方为0.862，模型拟合显著性检验的P值为0.00307 < 0.05，说明模型拟合效果还是可以的，
再看到参数拟合结果，x2的参数p值为0.728 > 0.05，因此下面我们剔除x2后再做一次拟合
'''

X1 = X[:,[0,1,3]]

# 用最小二乘法拟合回归方程
linear_model1 = sm.OLS(endog=y,exog=X1)
ols_result1 = linear_model1.fit()

# 输出拟合结果
print(ols_result1.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.901
Model:                            OLS   Adj. R-squared:                  0.862
Method:                 Least Squares   F-statistic:                     22.78
Date:                Sun, 22 May 2022   Prob (F-statistic):            0.00307
Time:                        17:18:21   Log-Likelihood:                -58.036
No. Observations:                   8   AIC:                             122.1
Df Residuals:                       5   BIC:                             122.3
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const        865.7005    447.091      1.936      0.111    -283.583    2014.984
x1           397.5845     60.183      6.606      0.001     242.879     552.290
x2          1120.7729    323.747      3.462      0.018     288.554    1952.992
==============================================================================
Omnibus:                        4.593   Durbin-Watson:                   1.536
Prob(Omnibus):                  0.101   Jarque-Bera (JB):                1.483
Skew:                           1.049   Prob(JB):                        0.477
Kurtosis:                       3.219   Cond. No.                         20.8
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

可以看到，模型整体的拟合是没有变的，在显著性水平为0.05的条件下，x1和x2的参数显著性通过检验，但常数项的显著性不明显，因此下面我们剔除常数项再做一次拟合
'''

X2 = X[:,[1,3]]

# 用最小二乘法拟合回归方程
linear_model2 = sm.OLS(endog=y,exog=X2)
ols_result2 = linear_model2.fit()

# 输出拟合结果
print(ols_result2.summary())
'''
                                 OLS Regression Results                                
=======================================================================================
Dep. Variable:                      y   R-squared (uncentered):                   0.986
Model:                            OLS   Adj. R-squared (uncentered):              0.981
Method:                 Least Squares   F-statistic:                              210.6
Date:                Sun, 22 May 2022   Prob (F-statistic):                    2.77e-06
Time:                        17:22:47   Log-Likelihood:                         -60.274
No. Observations:                   8   AIC:                                      124.5
Df Residuals:                       6   BIC:                                      124.7
Df Model:                           2                                                  
Covariance Type:            nonrobust                                                  
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
x1           499.5470     35.188     14.197      0.000     413.446     585.648
x2          1502.1518    310.270      4.841      0.003     742.948    2261.356
==============================================================================
Omnibus:                        0.202   Durbin-Watson:                   1.923
Prob(Omnibus):                  0.904   Jarque-Bera (JB):                0.253
Skew:                          -0.258   Prob(JB):                        0.881
Kurtosis:                       2.297   Cond. No.                         10.5
==============================================================================
Notes:
[1] R2 is computed without centering (uncentered) since the model does not contain a constant.
[2] Standard Errors assume that the covariance matrix of the errors is correctly specified.

可以看大，模型显著性和参数显著性都得到了提升，最终我们模型估计结果为：y = 499.5470*x1 + 1502.1518*x2
'''

2.1.2 自变量之间有交互作用的回归分析案例

某牙膏制造企业要求销售部门根据市场调查，找出公司生产的牙膏销售量与销售价格、广告投入等之间的关系，从而预测出在不同价格和广告费用下的销售量。为此，销售部收集了过去30个销售周期公司生产的牙膏的销售量、销售价格、广告费用，以及同期其他厂家生产的同类牙膏的平均销售价格，数据如下：

试根据这些数据建立一个数学模型，分析牙膏销售量与其他因素的关系，为制定价格策略和广告投入策略提供数据依据

在购买同类产品的牙膏时，顾客会在意不同品牌之间的价格差异，而不是价格本身。设牙膏销售量为y，价格差为x1，广告费为x2，其他厂家平均价格为x3，公司销售价格为x4，x1=x3-x4。

回归分析过程如下：

• （1）计算相关系数矩阵，观察自变量y与各因变量的相关程度，以及各因变量之间是否存在相关性；
• （2）绘制自变量与因变量的散点图，通过散点图初步判断模型是线性还是非线性；
• （3）建立模型，并对模型及其参数进行评估；
• （4）模型改进。在初始假定中，价格差和广告费对牙膏销售量的影响是相互独立的，而据直觉和经验可以猜想，价格差与广告费的交互作用也会对牙膏销售量有影响。
• （5）模型最终确定。

python代码如下：

import pandas as pd
import numpy as np
import statsmodels.api as sm

from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']    # 解决中文显示问题
plt.rcParams['axes.unicode_minus'] = False      # 解决坐标轴负数的负号显示问题

data = pd.read_excel(r'G:\牙膏销售数据表.xlsx')
print(data.info())

# 1、计算相关系数矩阵
print(data[['公司销售价格(元)','其他厂家平均价格(元)','广告费用(百万元)','价格差(元)','销售量(百万支)']].corr())
'''
             公司销售价格(元)  其他厂家平均价格(元)  广告费用(百万元)    价格差(元)  销售量(百万支)
公司销售价格(元)     1.000000     0.078367  -0.468793 -0.322067 -0.469220
其他厂家平均价格(元)   0.078367     1.000000   0.604540  0.918566  0.740948
广告费用(百万元)    -0.468793     0.604540   1.000000  0.759964  0.875954
价格差(元)       -0.322067     0.918566   0.759964  1.000000  0.889672
销售量(百万支)     -0.469220     0.740948   0.875954  0.889672  1.000000

可以看到，销售量与价格差、广告费用的相关系数分别是89%、88%，相关度很高，另外价格差和广告费用的相关度也不低（76%左右）
'''

# 2、绘制自变量与因变量的散点图
fig1 = plt.figure()
plt.scatter(x=data['广告费用(百万元)'].values,y=data['销售量(百万支)'].values,marker='*')
plt.xlabel('广告费用(百万元)')
plt.ylabel('销售量(百万支)')
fig2= plt.figure()
plt.scatter(x=data['价格差(元)'].values,y=data['销售量(百万支)'].values,marker='o')
plt.xlabel('价格差(元)')
plt.ylabel('销售量(百万支)')
plt.show()
'''
可以看到，价格差对销售量的走势是近似线性的，而广告费对销售量的走势有一点弯曲，
因此在假设价格差和广告费对牙膏销售量的影响是相互独立的前提下，准对价格差和销售量建立线性回归方程，准对广告费和销售量建立非线性回归方程（二次项）
'''

# 3、建立模型，并对模型及其参数进行估计
# 3.1、价格差和销售量的线性回归方程估计
X1 = sm.add_constant(data=data['价格差(元)'].values)    # X1由1列常数项和价格差组成
y = data['销售量(百万支)']

# 用最小二乘法估计方程
model1 = sm.OLS(endog=y,exog=X1)
result1 = model1.fit()
print(result1.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:               销售量(百万支)   R-squared:                       0.792
Model:                            OLS   Adj. R-squared:                  0.784
Method:                 Least Squares   F-statistic:                     106.3
Date:                Fri, 27 May 2022   Prob (F-statistic):           4.88e-11
Time:                        22:26:27   Log-Likelihood:                -7.0261
No. Observations:                  30   AIC:                             18.05
Df Residuals:                      28   BIC:                             20.85
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          7.8141      0.080     97.818      0.000       7.650       7.978
x1             2.6652      0.258     10.310      0.000       2.136       3.195
==============================================================================
Omnibus:                        5.481   Durbin-Watson:                   2.414
Prob(Omnibus):                  0.065   Jarque-Bera (JB):                4.092
Skew:                           0.883   Prob(JB):                        0.129
Kurtosis:                       3.391   Cond. No.                         4.69
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

可以看到，模型估计的拟合优度为0.792，在假设显著性水平α=0.05的条件下，模型显著性检验和参数显著性检验均通过，得到的回归方程为：y = 7.8141 + 2.6652*x1
'''

# 3.2、广告费和销售量的非线性回归方程估计
# 多项式处理
from sklearn.preprocessing import PolynomialFeatures
X2 = PolynomialFeatures(degree=2).fit_transform(X=data[['广告费用(百万元)']])

# 用最小二乘法估计方程
model2 = sm.OLS(endog=y,exog=X2)
result2 = model2.fit()
print(result2.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:               销售量(百万支)   R-squared:                       0.838
Model:                            OLS   Adj. R-squared:                  0.826
Method:                 Least Squares   F-statistic:                     69.81
Date:                Fri, 27 May 2022   Prob (F-statistic):           2.14e-11
Time:                        22:45:32   Log-Likelihood:                -3.2455
No. Observations:                  30   AIC:                             12.49
Df Residuals:                      27   BIC:                             16.69
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         25.1091      6.863      3.659      0.001      11.028      39.190
x1            -6.5589      2.217     -2.958      0.006     -11.109      -2.009
x2             0.6101      0.178      3.432      0.002       0.245       0.975
==============================================================================
Omnibus:                        0.063   Durbin-Watson:                   1.523
Prob(Omnibus):                  0.969   Jarque-Bera (JB):                0.224
Skew:                          -0.090   Prob(JB):                        0.894
Kurtosis:                       2.617   Cond. No.                     5.98e+03
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 5.98e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

可以看到，模型估计的拟合优度为0.838，在假设显著性水平α=0.05的条件下，模型显著性检验和参数显著性检验均通过，得到的回归方程为：y = 25.1091 - 6.5589*x2 + 0.6101*(x2)**2
'''

# 4、模型改进：加入交互项和广告费的二次项
data['价格差*广告费用'] = data['价格差(元)'] * data['广告费用(百万元)']
X3 = np.column_stack(tup=(X1,X2[:,[1,2]],data['价格差*广告费用'].values))

# 用最小二乘法估计方程
model3 = sm.OLS(endog=y,exog=X3)
result3 = model3.fit()
print(result3.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:               销售量(百万支)   R-squared:                       0.921
Model:                            OLS   Adj. R-squared:                  0.908
Method:                 Least Squares   F-statistic:                     72.78
Date:                Fri, 27 May 2022   Prob (F-statistic):           2.11e-13
Time:                        23:11:51   Log-Likelihood:                 7.5137
No. Observations:                  30   AIC:                            -5.027
Df Residuals:                      25   BIC:                             1.979
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         29.1133      7.483      3.890      0.001      13.701      44.525
x1            11.1342      4.446      2.504      0.019       1.978      20.291
x2            -7.6080      2.469     -3.081      0.005     -12.693      -2.523
x3             0.6712      0.203      3.312      0.003       0.254       1.089
x4            -1.4777      0.667     -2.215      0.036      -2.852      -0.104
==============================================================================
Omnibus:                        0.242   Durbin-Watson:                   1.512
Prob(Omnibus):                  0.886   Jarque-Bera (JB):                0.148
Skew:                          -0.153   Prob(JB):                        0.929
Kurtosis:                       2.843   Cond. No.                     9.81e+03
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 9.81e+03. This might indicate that there are
strong multicollinearity or other numerical problems.

可以看到，模型估计的拟合优度提升至0.921，在假设显著性水平α=0.05的条件下，模型显著性检验和参数显著性检验均通过，最终确定回归方程为：
y = 29.1133 + 11.1342*x1 - 7.6080*x2 + 0.6712*(x2)**2 - 1.4777*x1*x2
'''

2.1.3 非线性回归分析——预测第三产业国内生产总值案例

某省为了研究第三产业在本省宏观经济发展中的运行情况，对影响第三产业的资本要素、劳动力要素及科技进步要素这三项主要因素进行了统计分析，并运用道格拉斯生产函数建立了基本数学模型。原始数据如下：

现需要预测当资本投入为11738245、劳动力投入为987.37时，第三产业国内生产总值是多少。

道格拉斯生产函数对应数学模型为$Y=A{K}^{\alpha }{L}^{\beta }$，其中：

• Y是第三产业国内生产总值
• K是资金投入
• L是劳动力投入
• A是科技进步水平
• $\alpha$是资本弹性系数
• $\beta$是劳动弹性系数

上面的非线性模型我们可以通过对等号两边取对数，使之转换为多元线性模型，即$LnY=LnA+\alpha LnK+\beta LnL$，令$LnY=y，LnA=c，LnK=x_1，LnL=x_2$，于是得到$y=c+\alpha x_1+\beta x_2$

下面我们用python实现上述的建模过程：

import pandas as pd
import numpy as np
import statsmodels.api as sm

data = pd.read_excel(r'G:\第三产业国内生产总值数据表.xlsx')

# 线性化处理
data1 = data[['第三产业国内生产总值', '资本投入', '从业人员']].apply(lambda x:np.log(x))

# 确定自变量和因变量
X = sm.add_constant(data=data1[['资本投入', '从业人员']].values)  # 添加常数列
y = data1['第三产业国内生产总值'].values

# 用最小二乘法估计方程
model = sm.OLS(endog=y,exog=X)
result = model.fit()
print(result.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.988
Model:                            OLS   Adj. R-squared:                  0.985
Method:                 Least Squares   F-statistic:                     338.6
Date:                Tue, 31 May 2022   Prob (F-statistic):           1.86e-08
Time:                        23:20:30   Log-Likelihood:                 15.034
No. Observations:                  11   AIC:                            -24.07
Df Residuals:                       8   BIC:                            -22.87
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         -4.4745      1.229     -3.641      0.007      -7.308      -1.641
x1             0.4797      0.095      5.060      0.001       0.261       0.698
x2             0.6947      0.393      1.768      0.115      -0.212       1.601
==============================================================================
Omnibus:                        2.430   Durbin-Watson:                   1.134
Prob(Omnibus):                  0.297   Jarque-Bera (JB):                0.971
Skew:                           0.197   Prob(JB):                        0.615
Kurtosis:                       1.598   Cond. No.                         967.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
'''

根据返回的回归分析报告，我们可以看到：

• 样本回归方程为：$LnY=?4.4745+0.4797LnK+0.6947LnL$。
• 模型拟合情况：调整R方为0.985，接近于1，说明回归线对数据的拟合程度很高，拟合优度检验通过。
• 回归方程的总体显著性检验（F检验）：F统计量为338.6，对应的p值为$1.86\times 10^{?8}<0.05$，拒绝原假设，因此我们认为在0.05的显著性水平下（置信度为95%），模型的线性显著性通过检验。
• 回归系数的显著性检验（t检验）：常数项系数的t统计量对应的p值为0.007<0.05，LnK项系数的t统计量对应的p值为0.001<0.05，LnL项系数的t统计量对应的p值为0.115>0.05，在0.05的显著性水平下，我们可以拒绝常数项系数和LnK项系数为0的原假设，而没有充分的理由拒绝LnL项系数为0的原假设，在牺牲一定的置信度情况下，我们将显著性水平调整为0.12，从而各项系数均能通过显著性检验。

解释模型的经济含义：从回归方程可以看到，1992——2002年某省第三产业产出的资本投入弹性和劳动投入弹性分别是0.4797和0.6947。换而言之，在研究期间，保持劳动投入不变，每增加1个单位的资本投入，第三产业国内生产总值增加0.4797；同理，保持资本投入不变，每增加1个单位的劳动投入，第三产业国内生产总值增加0.6947。把两个产出弹性相加得到1.1744，就是规模报酬参数的取值，而该值大于1，说明某省第三产业经历了轻微的规模报酬递增。

模型预测：已知资本投入为11738245，劳动力投入为987.37，将其代入回归方程得到：$LnY=?4.4745+0.4797Ln(11738245)+0.6947Ln(987.37)=8.124218321936784$，即LnY=8.124218321936784，两边取自然指数e，得到Y=e^(8.124218321936784)=3375.2285581155024，因此当资本投入为11738245、劳动力投入为987.37时，第三产业国内生产总值是3375.23。

2.2 时间序列分析与市场预测

2.2.1 简单指数平滑（SES，Simple Exponential Smoothing）预测时间序列

该预测方法是适合于平稳序列（没有趋势和季节变动的序列），对过去的观测值加权平均进行预测的一种方法。观测值时间越远，其权数也跟着呈现指数的下降，因而称为指数平滑。

t+1的预测值是t期观测值与t期平滑值$S_t$的线性组合，其预测模型为：$F_{t+1}=\alpha Y_t+(1?\alpha )S_t$，其中

• $F_{t+1}$为t+1期的简单指数平滑预测值
• $Y_t$为第t期的实际观测值
• $S_t$为第t期的预测值，可设$S_0=Y_0$
• $\alpha$为平滑系数，$\alpha <1$

因此模型又可以表示为${\hat{y}}_{t+1|t}=\alpha y_t+\alpha (1?\alpha )y_{t?1}+\alpha (1?\alpha {)}^2{y}_{t?2}+?$

简单指数平滑预测系数α的选择：

• （1）不同的α会对预测结果产生不同的影响，当时间序列有较大的随机波动，宜选较小的α；注重于近期的实际值时，宜选较大的α
• （2）选择α时，还应考虑预测误差。预测误差可以用均方误差来衡量。
• （3）确定α时，可选择几个进行预测，然后找出预测误差最小的作为最后的值。
• （4）有两种极端情况：
  • $\alpha =0$时，所有未来值的预测等于历史数据的平均值，称为平均值法
  • $\alpha =1$时，所有未来值的预测设置为最后一次观测的值，统计中称为朴素方法

python代码实现案例：

import pandas as pd
from statsmodels.tsa.holtwinters import SimpleExpSmoothing
from matplotlib import pyplot as plt

dt = pd.DataFrame(data={
    'year':list(range(1990,2006)),
    '棉花产量(万吨)':[450.77,567.50,450.84,373.93,434.10,476.75,420.33,460.27,450.10,382.88,441.73,532.35,491.62,485.97,632.35,571.42]
})

fig0 = plt.figure()
plt.plot(dt['棉花产量(万吨)'].values)
plt.show()
'''
可以看到，折线图的走势没有明显的趋势和季节性变动
'''

# 下面我们分别用α=0.2，α=0.6，允许statsmodels自动查找最优α这三种方式进行拟合预测
fit1 = SimpleExpSmoothing(endog=dt['棉花产量(万吨)'].values).fit(smoothing_level=0.2,optimized=False)  # 用α=0.2拟合
fit2 = SimpleExpSmoothing(endog=dt['棉花产量(万吨)'].values).fit(smoothing_level=0.6,optimized=False)  # 用α=0.6拟合
fit3 = SimpleExpSmoothing(endog=dt['棉花产量(万吨)'].values).fit()  # 不做设置，自动查找最优α

# 绘制折线图，查看拟合和未来3年的预测值位置
fig = plt.figure()
line1 = plt.plot(list(fit1.fittedvalues)+list(fit1.forecast(steps=3)),c='g',marker='o')
line2 = plt.plot(list(fit2.fittedvalues)+list(fit2.forecast(steps=3)),c='r',marker='o')
line3 = plt.plot(list(fit3.fittedvalues)+list(fit3.forecast(steps=3)),c='b',marker='o')
line4 = plt.plot(dt['棉花产量(万吨)'].values,c='y',marker='*')
plt.legend(labels = ['alpha=0.2','alpha=0.6','auto','data'],loc='best')    # 添加图例
plt.show()

2.2.2 用Holt指数平滑模型预测2006年的人均GDP

Holt指数平滑(二次指数平滑)预测模型，一般简称为Holt模型，适合于含有趋势成分（或有一定的周期成分）序列的预测。

Holt模型有两个平滑系数和三个方程：

• 1、平滑值方程：$S_t=\alpha Y_t+(1?\alpha )(S_{t?1}+T_{t?1})$，其中$\alpha$是平滑参数（$0<\alpha <1$），$S_t$是t期的指数平滑值，$S_{t?1}$是t-1期的指数平滑值，$T_{t?1}$是t-1期的趋势值。该方程实际上是对t期平滑值$S_t$的修正，它把上一期的趋势值$T_{t?1}$加到$S_{t?1}$上，这样可以消除因趋势而产生的滞后，使其尽可能接近实际观测值$Y_t$。
• 2、趋势项更新方程：$T_t=\gamma (S_t?S_{t?1})+(1?\gamma )T_{t?1}$，其中$\gamma$是平滑系数（$0<\gamma <1$），$T_t$是t期趋势的平滑值。该方程实际上是对趋势的修正，$T_t$被表示成相邻两项平滑值之差（$S_t?S_{t?1}$），如果序列存在趋势，则新的观测值总是高于（上升趋势）或低于（下降趋势）前一期数值，同时由于随机波动的影响，需要用$\gamma$平滑（$S_t?S_{t?1}$）的趋势，然后再将平滑的结果加到前一期趋势的估计值$T_{t?1}$与$1?\gamma$的乘积上。
• K期预测值方程：$F_{t+k}=S_t+k\times T_t$，其中k是用于预测的时期数，当k=1时，t+1期的预测值就是t期的平滑值$S_t$加上t期的修正趋势值$T_t$。

Holt模型中初始值的确定：

• （1）由于在开始计算时，还没第1个时期的平滑值$S_1$和修正趋势值$T_1$，通常设$S_1=Y_1，T_1=Y_2?Y_1$；
• （2）平滑系数$\alpha 和\gamma$可根据实际情况确定，可用均方误差来衡量。

案例

1990年到2005年的人均GDP数据如下：

可以看到人均GDP逐年增加，存在增长趋势，下面我们利用Holt模型预测2006年的人均GDP：

• （1）由于我们要预测的是下一期（即2006年）的人均GDP，所以k=1，于是K期预测值方程转换为$F_{t+1}=S_t+T_t$；
• （2）在开始计算时，我们令$S_1=Y_1，T_1=Y_2?Y_1$，于是对应地可以得到$S_{1990}=Y_{1990}，T_{1990}=Y_{1991}?Y_{1990}$；
• （3）确定好初始值之后，我们来确定平滑系数$\alpha$和$\gamma$，一般我们认为过去历史值对现在的影响比较弱小，下一期的变化走势主要由本期决定，因此可以设置平滑系数为0.7；
• （4）确定好各项系数后，我们就可以利用预测方程$F_{t+1}=S_t+T_t=0.7\times Y_t+0.3\times (S_{t?1}+T_{t?1})+0.7\times (S_t?S_{t?1})+0.3\times T_{t?1}$进行迭代：
  • （1）$F_{1991}=S_{1990}+T_{1990}=Y_{1990}+Y_{1991}?Y_{1990}=Y_{1991}=1892.76，S_{1991}=0.7\times Y_{1991}+0.3\times (S_{1990}+T_{1990})=Y_{1991}=1892.76，T_{1991}=0.7\times (S_{1991}?S_{1990})+0.3\times T_{1990}=T_{1990}=Y_{1991}?Y_{1990}=248.29$
  • （2）$F_{1992}=S_{1991}+T_{1991}=1892.76+248.29=2141.05，S_{1992}=0.7\times Y_{1992}+0.3\times (S_{1991}+T_{1991})=2260.078，T_{1992}=0.7\times (S_{1992}?S_{1991})+0.3\times T_{1991}=331.6096$
  • （3）以此类推，最后得到$F_{2006}=S_{2005}+T_{2005}=15589.21213$

最终得到数据如下：

对于长期预测，使用Holt方法的预测在未来会无限期地增加或减少，在这种情况下，我们使用具有阻尼参数$?(0<?<1)$的阻尼趋势方法来防止预测 “失控”，因此我们对上面三个方程进行优化：

• 平滑值方程：$S_t=\alpha Y_t+(1?\alpha )(S_{t?1}+?T_{t?1})$
• 趋势项方程：$T_t=\gamma (S_t?S_{t?1})+(1?\gamma )?T_{t?1}$
• 预测值方程：$F_{t+k}=S_t+(?+{?}^2+?+{?}^k)T_t$

python代码实现案例：

import pandas as pd
from statsmodels.tsa.holtwinters import Holt
from matplotlib import pyplot as plt

dt = pd.DataFrame(data={
    'year':list(range(1990,2006)),
    '人均GDP':[1644.47,1892.76,2311.09,2998.36,4044.00 ,5045.73,5845.89,6420.18,6796.03,7158.50,7857.68,8621.71,9398.05,10541.97,12335.58,14040.00]
})

fig0 = plt.figure()
plt.plot(dt['人均GDP'].values)
plt.show()
'''
可以看到，折线图的走势有明显的上升趋势
'''

# 下面我们尝试Holt方法的三种变体，对比拟合效果
fit1 = Holt(endog=dt['人均GDP'].values).fit(smoothing_level=0.8,smoothing_trend=0.2,optimized=False)  # 设置alpha=0.8，gamma=0.2
fit2 = Holt(endog=dt['人均GDP'].values,exponential=True).fit(smoothing_level=0.8,smoothing_trend=0.2,optimized=False)  # 使用指数模型，而不是Holt的加法模型（默认）
fit3 = Holt(endog=dt['人均GDP'].values,damped_trend=True).fit(smoothing_level=0.8,smoothing_trend=0.2)  # 使用阻尼版本的Holt附加模型，但允许优化阻尼参数phi，同时固定alpha=0.8，gamma=0.2的值

# 绘制折线图，查看拟合和未来3年的预测值位置
fig = plt.figure()
line1 = plt.plot(list(fit1.fittedvalues)+list(fit1.forecast(steps=3)),c='g',marker='.')
line2 = plt.plot(list(fit2.fittedvalues)+list(fit2.forecast(steps=3)),c='r',marker='.')
line3 = plt.plot(list(fit3.fittedvalues)+list(fit3.forecast(steps=3)),c='b',marker='.')
line4 = plt.plot(dt['人均GDP'].values,c='y',marker='^')
plt.legend(labels = ["Holt's linear trend",'Exponential trend','Additive damped trend','data'],loc='best')    # 添加图例
plt.show()

2.2.3 指数曲线预测

当时间序列以几何级递增或递减时，适合用指数曲线对样本进行拟合，其模型的一般形式为${\hat{Y}}_t=b_0\exp (b_{1}t)=b_0{e}^{b_{1}t}$，式中$b_0、b_1$为待定系数，exp表示自然对数$\ln$的反函数，e=2.71828182845904。

对上面指数曲线模型线性化处理，即等号两边取自然对数得到：$\ln {\hat{Y}}_t=\ln b_0+b_{1}t$。将模型线性化之后，我们就可以使用最小二乘法来估计参数。

2.2.4 多项式曲线预测

有些现象的变化形态比较复杂，它们不是按照某种固定的形态变化，而是有升有降，在变化过程中可能有几个拐点，这时就需要拟合多项式函数：

• 当只有一个拐点时，可以拟合二阶曲线，即抛物线；
• 当有两个拐点时，需要拟合三阶曲线；
• 当有k-1个拐点时，需要拟合k阶曲线。

k阶曲线函数的一般形式为：$\widehat{Y_t}=b_0+b_{1}t+b_2{t}^2+b_k{t}^k$。

将函数线性化处理：令$t=x_1，t^2=x_2，?，t^k=x_k$。经过处理后，函数变成多元线性回归方程，于是我们就可以使用最小二乘法来估计$b_0，b_1，b_2，?，b_k$。

2.2.5 用Winters指数平滑模型预测啤酒的销售量

简单指数平滑模型适合于对平稳序列（没有趋势和季节成分）的预测；Holt指数平滑模型适合于含有趋势成分但不含季节成分序列的预测。

如果时间序列中既含有趋势成分又含有季节成分，则可以使用Winters指数平滑模型进行预测，该模型要求数据是按季度或月份收集的，而且至少需要4年（4个季节周期长度）以上的数据。

Winters指数平滑(三次指数平滑)模型包括三个平滑参数$\alpha 、\gamma 、\delta$（取值均在0和1之间）和四个方程：

• （1）平滑值：$S_t=\alpha \frac{Y_t}{I_{t?L}}+(1?\alpha )(S_{t?1}+T_{t?1})$，其中L为季节周期的长度，若为季节数据，则L=4，若为月份数据，则L=12；I为季节调整因子，$\frac{Y_t}{I_{t?L}}$表示t期观测值剔除季节调整因子$I_{t?L}$来消除季节变动。
• （2）趋势项更新：$T_t=\gamma (S_t?S_{t?1})+(1?\gamma )T_{t?1}$，用参数$\gamma$加权趋势增值$(S_t?S_{t?1})$，用$(1?\gamma )$加权前期趋势值，以此来对趋势值$T_{t?L}$进行修正。
• （3）季节项更新：$I_t=\delta \frac{Y_t}{S_t}(1?\delta )I_{t?L}$，其中$\frac{Y_t}{S_t}$是根据季节变动来调整实际值，用参数$\delta$加权这一调整值，用$(1?\delta )$加权前一个季度数据$I_{t?L}$，其结果就是t期的季节调整因子。
• （4）K期预测：$F_{t+k}=(S_t+k{T}_t)I_{t?L+k}$。

在进行预测时，模型中的参数即$\alpha 、\gamma 、\delta$难以确定时，可采用spss软件自动搜寻的方式确定，而初始平滑值和趋势值也可以采用spss软件自动的方式确定。

Winters指数平滑法有两种变体：

• 加法形式：整个序列的季节变化基本保持不变。
• 乘法形式：季节变化与系列水平成比例变化。

python代码实现案例：

import pandas as pd
from statsmodels.tsa.holtwinters import ExponentialSmoothing
from matplotlib import pyplot as plt

dt = pd.DataFrame(data={
    'year_season':['2005_Q1','2005_Q2','2005_Q3','2005_Q4','2006_Q1','2006_Q2','2006_Q3','2006_Q4','2007_Q1','2007_Q2','2007_Q3','2007_Q4','2008_Q1','2008_Q2','2008_Q3','2008_Q4','2009_Q1','2009_Q2','2009_Q3','2009_Q4','2010_Q1','2010_Q2','2010_Q3','2010_Q4'],
    '销售量':[25,32,37,26,30,38,42,30,29,39,50,35,30,39,51,37,29,42,55,38,31,43,54,41]
})

fig0 = plt.figure()
plt.plot(dt['销售量'].values)
plt.show()
'''
可以看到，折线图的走势有明显的上升趋势和季节性变动
'''

# 下面我们尝试完整的Holt-Winters方法，包括趋势组件和季节性组件，对比拟合效果
fit1 = ExponentialSmoothing(endog=dt['销售量'].values,trend='add',seasonal='add',seasonal_periods=4).fit()  # 使用加法趋势，周期season_length = 4的加性季节
fit2 = ExponentialSmoothing(endog=dt['销售量'].values,trend='add',seasonal='mul',seasonal_periods=4).fit()  # 使用加法趋势，周期season_length = 4的乘法季节
fit3 = ExponentialSmoothing(endog=dt['销售量'].values,trend='add',seasonal='add',damped_trend=True,seasonal_periods=4).fit()  # 使用加性阻尼趋势，周期season_length = 4的加性季节
fit4 = ExponentialSmoothing(endog=dt['销售量'].values,trend='add',seasonal='mul',damped_trend=True,seasonal_periods=4).fit()  # 使用加性阻尼趋势，周期season_length = 4的加性季节

# 绘制折线图，查看拟合和未来3年的预测值位置
fig = plt.figure()
line1 = plt.plot(list(fit1.fittedvalues)+list(fit1.forecast(steps=4)),c='g',marker='*')
line2 = plt.plot(list(fit2.fittedvalues)+list(fit2.forecast(steps=4)),c='r',marker='*')
line3 = plt.plot(list(fit3.fittedvalues)+list(fit3.forecast(steps=4)),c='b',marker='*')
line4 = plt.plot(list(fit3.fittedvalues)+list(fit3.forecast(steps=4)),c='b',marker='*')
line5 = plt.plot(dt['销售量'].values,c='y',marker='^')
plt.legend(labels = ["aa",'am','aa damped','am damped','data'],loc='best')    # 添加图例
plt.show()

2.2.6 分解预测法预测啤酒的销售量

分解预测法的步骤：

• 1、计算季节指数：以其平均数等于100%为条件而构成的反映季节变动的值，表示某一月份或季度的数值占全年平均数值的大小。如果现象的发展没有季节变动，则各期的季节指数应等于100%。季节变动的程度是根据各季节指数与其平均数（100%）的偏差程度来测定。季节指数计算步骤如下：
  • （1）计算移动平均值（季节数据采用4项移动平均，月份数据采用12项移动平均），并将其结果进行 “中心化” 处理（即将移动平均的结果再进行一次二项移动平均）。
  • （2）计算移动平均的比值（也称为季节比率）：将序列的各项观测值除以相应的中心化移动平均值。然后再计算出各比值的季度（或月份）平均值，即季节指数。
  • （3）季节指数调整：各季节指数的平均数应等于1或100%，若根据上面计算的季节比率的平均值不等于1时，则需要进行调整，具体方法是将上面计算得出的每个季节比率除以它们的总平均值。
• 2、确定并分离季节成分：分离季节成分，即将原时间序列除以相应的季节指数。季节因素分离后的序列反映了在没有季节因素影响的情况下时间序列的变化形态：$\frac{Y}{S}=\frac{T\times S\times I}{S}=T\times I$。
• 3、建立线性预测模型进行预测：
  • （1）根据分离季节性因素的序列确认线性趋势方程。
  • （2）根据趋势方程进行预测，该预测值不含季节性因素，即为没有季节因素影响情况下的预测值。
• 4、计算出最后的预测值：将回归预测值乘以相应的季节指数。

案例：

（1）计算移动平均和季节比率

（2）季节指数调整

（3）分离季节成分，建立线性预测模型进行预测

线性趋势模型建立：

import pandas as pd
import statsmodels.api as sm

dt = pd.DataFrame(data={
    'time':list(range(1,25)),
    '分离季节成分':[31.48190746,30.62685489,29.12459193,29.19202024,37.77828896,36.36939019,33.06034759,33.68310027,36.51901266,37.3264794,39.35755666,39.29695032,37.77828896,37.3264794,40.14470779,41.54249034,36.51901266,40.19774705,43.29331232,42.66526034,39.03756526,41.15483626,42.50616119,46.03357037]
})

X = sm.add_constant(data=dt['time'].values)  # 添加常数列
y = dt['分离季节成分'].values

# 用最小二乘法估计分离季节成分后的线性趋势方程
model = sm.OLS(endog=y,exog=X)
result = model.fit()
print(result.summary())
'''
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.767
Model:                            OLS   Adj. R-squared:                  0.756
Method:                 Least Squares   F-statistic:                     72.22
Date:                Sun, 19 Jun 2022   Prob (F-statistic):           2.13e-08
Time:                        12:08:25   Log-Likelihood:                -52.328
No. Observations:                  24   AIC:                             108.7
Df Residuals:                      22   BIC:                             111.0
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         30.5780      0.942     32.449      0.000      28.624      32.532
x1             0.5605      0.066      8.498      0.000       0.424       0.697
==============================================================================
Omnibus:                        0.733   Durbin-Watson:                   1.698
Prob(Omnibus):                  0.693   Jarque-Bera (JB):                0.695
Skew:                          -0.089   Prob(JB):                        0.707
Kurtosis:                       2.186   Cond. No.                         29.6
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
'''

从最小二乘法回归结果可以看到，模型拟合的调整R方不是很高，但模型F检验和参数的t检验在显著性水平$\alpha$为0.05的情况下，均通过检验，因此我们确定分离季节成分后的序列的线性趋势方程为$\widehat{Y_t}=30.5780+0.5605t$。

（4）根据线性趋势方程进行预测，并计算最后的预测值

tips：以上大部分知识点主要来自书籍《CDA数据分析实务》的总结，同时也借鉴了文章《如何使用Python构建指数平滑模型:Simple Exponential Smoothing, Holt, and Holt-Winters》http://t.zoukankan.com/harrylyx-p-11852149.html。