Python实现：时间序列趋势外推法应用-龚珀兹曲线拟合

龚珀兹曲线

下表数据为某跨国公司1989-2021年的年销售量数据，使用适合的模型预测该公司2022年的销售额，并得出理由。

部分数据如下表（具体数据从主页资源下载）：

1、读取数据： 

#读取数据
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

data = pd.read_excel('ch4综合分析.xlsx')
data.head()
#%%
# 创建散点图
#显示中文
from matplotlib import rcParams

# 设置中文字体
rcParams['font.sans-serif'] = ['SimHei']  # 用黑体显示中文
rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

plt.scatter(data.index, data['总额（yt）'], color='blue', label='原始数据')
plt.xlabel('时序（t）')
plt.ylabel('总额（yt）')
plt.legend()
plt.show()

 通过散点图得知数据符合龚珀兹曲线特征，为此采用龚珀兹曲线拟合。

2、原理

（1）取对数：

（2）分组法求解参数

或者

（3）Python求解

#选择合适的趋势外推法预测销售额
#拟合
# 对 '总额（yt）' 列取对数并写入表格
data['总额（yt）_log'] = np.log10(data['总额（yt）'])
# 显示前几行数据
print(data.head(33))
print('----------------------------------------------------')
#使用分组法估计参数
#分组
n = (2021-1989+1)/3
#读取前n个数据
data1 = data.iloc[:11]
data2 = data.iloc[11:22]
data3 = data.iloc[22:]
#分组求和
data1_sum = data1['总额（yt）_log'].sum()
data2_sum = data2['总额（yt）_log'].sum()
data3_sum = data3['总额（yt）_log'].sum()
#计算参数
b = ((data3_sum-data2_sum)/(data2_sum-data1_sum))**(1/11)
a = (data2_sum-data1_sum)*((b-1)/(b**11-1)**2)
k = (data1_sum-a*(b**11-1)/(b-1))*(1/11)

a1=10**a
k1=10**k
#输出参数
print('k1 = %f, a1 = %f, b = %f' % (k1,a1,b))

#输出拟合函数
print('拟合函数为：y = %f * %f ^ %f^t' % (k1,a1,b))

#计算第2022年的销售量
t1=33
d34 = k1 *a1 ** b**t1
print('2022年的销售量为：%d' % d34)
#计算se
# 计算预测值
data['预测值'] = k1 * a1 ** (b ** data.index)
# 计算残差
data['残差'] = data['总额（yt）'] - data['预测值']
# 计算残差的标准误差
se = np.sqrt(np.sum(data['残差'] ** 2) / (len(data) - 2))
# 输出标准误差
print('标准误差 (SE) = %f' % se)

（4）龚珀兹曲线方程

3、曲线拟合 

import matplotlib.pyplot as plt
# 计算拟合值
data['拟合值'] = k1 * a1 ** (b ** data.index)
# 创建散点图和拟合曲线
plt.scatter(data.index, data['总额（yt）'], color='blue', label='原始数据')
plt.plot(data.index, data['拟合值'], color='red', label='拟合曲线')
#添加预测点
plt.scatter(t1, d34, color='green', label='2022年预测点')
plt.xlabel('时序（t）')#时序从第0开始计时
plt.ylabel('总额（yt）')
plt.legend()
plt.show()