基于线性回归的商品销售额预测（假日效应、季节性商品处理）

(一)环境安装及配置

“holidays” 是一个Python包，用于处理节假日和工作日的计算和管理。它提供了一种方便的方式来识别不同国家或地区的公共假日，以及判断指定日期是否为工作日。这在进行日期相关的分析和建模时非常有用，例如在预测销售额、交通流量等方面。

!pip uninstall holidays -y
!pip install holidays==0.18

import datetime as dt
import numpy as np
import pandas as pd
from sklearn import linear_model
from sklearn import metrics
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import GroupKFold
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
pd.options.mode.chained_assignment = None  # default='warn'
train = pd.read_csv('../input/playground-series-s3e19/train.csv')
test = pd.read_csv('../input/playground-series-s3e19/test.csv')
test['num_sold'] = 300 

import matplotlib
matplotlib.rcParams.update({'font.size': 15})

导入所需的库：

datetime as dt：用于处理日期和时间。
numpy as np：用于数值计算。
pandas as pd：用于处理和操作结构化数据。
linear_model 和 metrics 子模块从 sklearn：用于线性模型和评估指标。
matplotlib.pyplot as plt：用于绘制图表和可视化。
seaborn as sns：用于统计数据可视化。
GroupKFold、StandardScaler 和 Pipeline 从 sklearn.model_selection 和 sklearn.preprocessing：用于交叉验证、特征标准化和构建处理管道。

禁用 pandas 的链式赋值警告：
代码禁用了 pandas 对于链式赋值操作的警告提示，以避免显示警告信息。

读取训练数据和测试数据：
使用 pd.read_csv 从 CSV 文件中读取训练数据和测试数据，并分别存储为名为 train 和 test 的 pandas DataFrame。

添加测试数据的列：
在测试数据 test 中添加一个名为 ‘num_sold’ 的列，所有行的值都设置为 300。这是虚拟的列，用于后续的数据合并操作。
更新 matplotlib 字体大小：
通过 matplotlib.rcParams.update 更新图表的字体大小为 15，以改善图表的可读性。

(二) 数据EDA

利用 2017 年至 2021 年的销售数据，我们的目标是构建一个预测模型，用于预测即将到来的 2022 年的“num_sold”。这种预测依赖于对四个要素的全面了解：日期、国家/地区、商店和不同的产品类别。因此，我们的第一步是系统地分析“num_sold”与这四个参数中的每一个之间的关系。

train.head(5)

# Creating variables for analysis
analysis = train.drop(columns='id')
# date column is separated for each element
analysis['date'] = pd.to_datetime(analysis['date'])
analysis['day'] = analysis['date'].dt.day
analysis['week'] = analysis['date'].dt.dayofweek
analysis['month'] = analysis['date'].dt.month
analysis['year'] = analysis['date'].dt.year
analysis['day_of_year'] = analysis['date'].dt.dayofyear
analysis['time_no'] = (
    analysis['date'] - dt.datetime(2017, 1, 1)) // dt.timedelta(days=1)
analysis.loc[analysis['date'] > dt.datetime(2020, 2, 29), 'time_no'] -= 1
date_columns = ['date', 'day', 'week', 'month', 'year', 'time_no']

1. 删除列： 从DataFrame train 中删除了“id”列，并将结果存储在名为analysis的新DataFrame中。

2. 转换为日期时间： 使用pd.to_datetime()将analysis DataFrame中的“date”列转换为pandas的日期时间格式。

3. 提取日期组件： 使用pandas的日期时间功能提取了日期的各个组件，例如日、星期、月、年、年中的天数以及自定义的“time_no”值。

   • day：从“date”列提取月份中的日期。
   • week：从“date”列提取星期几（0表示星期一，6表示星期日）。
   • month：从“date”列提取月份。
   • year：从“date”列提取年份。
   • day_of_year：从“date”列提取一年中的天数（1到365/366）。
   • time_no：创建一个自定义的时间值，计算自2017年1月1日以来的天数。对2020年2月29日之后的日期进行了调整。

4. 调整闰年的time_no： 对于2020年2月29日之后的日期，从“time_no”值中减去了1。是为了考虑闰年中多出的一天。

5. 日期列： 定义了一个名为date_columns的列表，其中包含与日期信息相关的列的名称。

uniques = {}
for column in analysis.columns:
    uniques[column] = analysis[column].unique().tolist()
print(uniques['country'])
print(uniques['store'])
print(uniques['product'])

遍历analysis DataFrame 的每个列，并创建一个字典 uniques，其中包含每列的唯一值列表。打印出 'country'、'store' 和 'product' 列的唯一值列表。

看一下数据中这些特定列的不同取值，从而检查数据质量、查找异常值。

2.1 分析存储数据

使用matplotlib和Seaborn库绘制了3个子图，进行销售数据的可视化。

fig, axs = plt.subplots(1, 3, figsize=(20, 5))

# First plot
grouped_data = analysis.groupby(['date', 'store'])['num_sold'].sum().reset_index()
for store in uniques['store']:
    store_data = grouped_data[grouped_data['store'] == store]
    axs[0].plot(store_data['date'], store_data['num_sold'], ".", label=store)
axs[0].legend()
axs[0].set_title('Aggregated Sales Over Time Per Store')

# Second plot
for store in uniques['store']:
    store_data = grouped_data[grouped_data['store'] == store]
    axs[1].plot(store_data['date'], store_data['num_sold'] /
             store_data['num_sold'].sum(), ".", label=store)
axs[1].legend()
axs[1].set_title('Aggregated Sales after Normalization Over Time Per Store')

# Third plot
for store in uniques['store']:
    store_data = grouped_data[grouped_data['store'] == store]['num_sold']
    sum_store=store_data.sum()
    axs[2] = sns.kdeplot(store_data/sum_store, label=store, ax=axs[2])
axs[2].legend()
axs[2].set_title('Normalized Sales Distribution Across Stores')

plt.tight_layout()
plt.show()

1. 第一个子图：聚合的销售趋势（按商店）

   • 该子图展示了不同商店的销售趋势随时间的变化。
   • 使用groupby函数将数据按日期和商店进行分组，计算每个商店在每天的销售总量。
   • 对于每个商店，通过在图中绘制点来表示每天的销售情况，点的颜色可能是不同商店之间的区分。
   • axs[0]表示第一个子图，legend()添加了图例，set_title()设置了图的标题

2. 第二个子图：归一化的聚合销售趋势（按商店）

   • 该子图显示了不同商店的归一化销售趋势随时间的变化。
   • 归一化处理是为了在不同商店之间进行比较，以便更好地理解它们的相对销售情况。
   • 归一化的销售量是每个商店的销售量除以该商店总销售量的结果。
   • 同样，使用点来表示每天的归一化销售情况，点的颜色区分不同商店。

3. 第三个子图：归一化销售分布（跨商店）

   • 该子图显示了归一化的销售量分布在不同商店之间的情况。
   • 使用Seaborn的kdeplot绘制了销售量的核密度估计，以展示销售量分布的形状。
   • 每个商店的核密度估计曲线叠加在一起，通过图例区分不同商店。

在检查这些商店时，似乎每个商店都显示相似的销售分布，只是通过标量常量来区分。如上所示，当按总销售额对每个商店的销售额进行归一化并比较分布时，会出现这种模式

2.2 分析国家数据

fig, axs = plt.subplots(1, 3, figsize=(20, 5))

# First plot
grouped_data_country = analysis.groupby(['date', 'country'])['num_sold'].sum().reset_index()


for country in uniques['country']:
    country_data = grouped_data_country[grouped_data_country['country'] == country]
    axs[0].plot(country_data['date'], country_data['num_sold'], ".", label=country)
axs[0].legend()
axs[0].set_title('Aggregated Sales Over Time Per Country')

# Second plot
for country in uniques['country']:
    country_data = grouped_data_country[grouped_data_country['country'] == country]
    axs[1].plot(country_data['date'], country_data['num_sold']/country_data['num_sold'].sum(), ".", label=country)
axs[1].legend()
axs[1].set_title('Aggregated Sales after Normalization Over Time Per Country')

# Third plot
for country in uniques['country']:
    country_data = grouped_data_country[grouped_data_country['country']== country]
    sum_country = country_data['num_sold'].sum()
    axs[2] = sns.kdeplot(country_data['num_sold']/sum_country, label=country)
axs[2].legend()
axs[2].set_title('Normalized Sales Distribution Across Countries')
plt.tight_layout()
plt.show()

1. 第一个子图：聚合的销售趋势（按国家）

   • 该子图显示了不同国家的销售趋势随时间的变化。
   • 使用groupby函数将数据按日期和国家进行分组，计算每个国家在每天的销售总量。
   • 对于每个国家，通过在图中绘制点来表示每天的销售情况，点的颜色可能是不同国家之间的区分。
   • axs[0]表示第一个子图，legend()添加了图例，set_title()设置了图的标题。

2. 第二个子图：归一化的聚合销售趋势（按国家）

   • 该子图显示了不同国家的归一化销售趋势随时间的变化。
   • 归一化处理是为了在不同国家之间进行比较。
   • 归一化的销售量是每个国家的销售量除以该国家总销售量的结果。
   • 同样，使用点来表示每天的归一化销售情况，点的颜色区分不同国家。

3. 第三个子图：归一化销售分布（跨国家）

   • 该子图显示了归一化的销售量分布在不同国家之间的情况。
   • 使用Seaborn的kdeplot绘制了销售量的核密度估计，以展示销售量分布的形状。
   • 每个国家的核密度估计曲线叠加在一起，通过图例区分不同国家。

tips：年份差异是由于不同国家的GDP差异造成的，也可以将GDP用作外部预测因子。

import requests
def get_gdp_per_capita(country, year):
    alpha3 = {'Argentina': 'ARG', 'Canada': 'CAN',
              'Estonia': 'EST', 'Japan': 'JPN', 'Spain': 'ESP'}
    url = "https://api.worldbank.org/v2/country/{0}/indicator/NY.GDP.PCAP.CD?date={1}&format=json".format(
        alpha3[country], year)
    response = requests.get(url).json()
    return response[1][0]['value']

gdp = []
for country in train.country.unique():
    row = []
    for year in range(2017, 2023):
        row.append(get_gdp_per_capita(country, year))
    gdp.append(row)
    
gdp = np.array(gdp)
gdp /= np.sum(gdp) #to renomralize the data

rel_gdp_df = pd.DataFrame(gdp, index=analysis.country.unique(), columns=range(2017, 2023))

df= analysis.groupby(['date', 'country'])[['num_sold']].sum().reset_index()
df['rel_gdp'] = df.apply(lambda s: rel_gdp_df.loc[s.country, s.date.year], axis=1)

fig, axs = plt.subplots(1, 2, figsize=(15, 5))

for country in uniques['country']:
    country_data = df[df['country'] == country]
    axs[0].plot(country_data['date'], country_data['num_sold']/country_data['rel_gdp'], ".", label=country)
axs[0].legend()
axs[0].set_title('Aggregated Sales divided by GDP Over Time Per Country')

# Third plot
for country in uniques['country']:
    country_data = df[df['country'] == country]
    axs[1] = sns.kdeplot(country_data['num_sold']/country_data['rel_gdp'], label=country)
axs[1].legend()
axs[1].set_title(
    'Normalized Sales Distribution  divided by GDP Across Countries')
plt.tight_layout()
plt.show()

每个国家在按GDP归一化后都有大致相同的分布。差异性可能来源于节日庆祝日期不同。
通过考虑假日以及将销售量与相对GDP进行比较，来深入分析销售数据：

# Generate holidays for each country and year
import holidays
for year in [2017, 2018, 2019, 2020,2021]:
    for country in uniques['country'][:]:
        for date, holiday_name in sorted(holidays.CountryHoliday(country, years=year).items()):
            df = df[~((df['date'].dt.date >= date) & (
                df['date'].dt.date < date+dt.timedelta(days=10)))]

fig, axs = plt.subplots(1, 2, figsize=(15, 5))

for country in uniques['country']:
    country_data = df[df['country'] == country]
    axs[0].plot(country_data['date'], country_data['num_sold'] /
                country_data['rel_gdp'],label=country)
axs[0].legend()
axs[0].set_title(f'Aggregated Sales divided by GDP Over Time Per Country\n (after removing holidays)')

# Third plot
for country in uniques['country']:
    country_data = df[df['country'] == country]
    axs[1] = sns.kdeplot(country_data['num_sold'] /
                         country_data['rel_gdp'], label=country)
axs[1].legend()
axs[1].set_title(
    'Normalized Sales Distribution  divided by GDP Across Countries')
plt.tight_layout()
plt.show()

1. 生成假日信息：

• 使用holidays库，为每个国家和年份（2017到2021）生成节假日信息。
• 通过嵌套循环，对于每个国家和每年，循环遍历假日并从数据中删除与假日日期相接近的10天范围内的数据行。

2. 绘制图形：

• 生成了两个子图 第一个子图绘制了每个国家的销售量除以相对GDP的趋势随时间的变化。
• 第二个子图使用Seaborn的kdeplot绘制了销售量除以相对GDP的核密度估计，以显示销售分布的形状。
  
  使用任何一个 apporach 来规范化数据，但第二个看起来更稳定。

2.3 产品数据分析

# First plot
fig, axs = plt.subplots(1, 2, figsize=(15, 5))

• 创建一个包含两个子图的图像，每个子图的大小为 (15, 5)。

grouped_data = analysis.groupby(['date', 'country','product'])['num_sold'].sum().reset_index()

• 使用 groupby 对原始数据进行分组，以日期、国家和产品作为分组依据，计算每个组的销售总量，并将结果重置索引，得到一个新的 DataFrame。

grouped_data['rel_gdp'] = grouped_data .apply(lambda s: rel_gdp_df.loc[s.country, s.date.year], axis=1)

• 为上述分组后的 DataFrame 添加一个相对 GDP 列。通过 apply 函数，根据每行的国家和年份获取相对 GDP 的值。

grouped_data['num_sold'] = grouped_data['num_sold'] / grouped_data['rel_gdp']

• 将销售量除以相对 GDP，以进行销售量的归一化。

grouped_data = grouped_data.groupby(['date','product'])['num_sold'].sum().reset_index()

• 再次使用 groupby 对数据进行分组，以日期和产品作为分组依据，计算每个日期中每个产品的销售总量。

for product in uniques['product']:
    product_data = grouped_data[grouped_data['product'] == product]
    axs[0].plot(product_data['date'], product_data['num_sold'],label=product)

• 在第一个子图中，对每个产品循环遍历。
• 对于每个产品，从 grouped_data 中选择匹配产品的数据。
• 使用 axs[0].plot 在第一个子图上绘制每个产品的销售趋势。

axs[0].set_title('Aggregated Sales Over Time Per Product')

• 为第一个子图设置标题。

grouped_data = analysis.groupby(['date', 'product'])['num_sold'].sum().reset_index()

• 使用 groupby 对原始数据进行分组，以日期和产品作为分组依据，计算每个日期中每个产品的销售总量，并将结果重置索引，得到一个新的 DataFrame。

for product in uniques['product']:
    product_data = grouped_data[grouped_data['product'] == product]
    axs[1].plot(product_data['date'], product_data['num_sold']/product_data['num_sold'].sum(),label=product)

• 在第二个子图中，对每个产品循环遍历。
• 对于每个产品，从 grouped_data 中选择匹配产品的数据。
• 使用 axs[1].plot 在第二个子图上绘制每个产品的归一化销售趋势。

axs[1].legend(bbox_to_anchor=(1.05, 1), loc='upper left')

• 为第二个子图添加图例，并通过 bbox_to_anchor 和 loc 参数调整图例的位置。

axs[1].set_title('Aggregated Sales after Normalization Over Time Per Product')

• 为第二个子图设置标题。

plt.tight_layout()
plt.show()

• 调整子图布局并显示整个图像。
  
  对于每种产品，我们观察销售趋势，以展示正弦和余弦成分，指示周期性模式。然而，这些周期性成分的振幅因产品而异。因此，我们将分别对待产品。

此外，还有与假日季节一致的明显销售高峰。

由于新冠病毒，2020年是不同的一年。

2.4 分析新冠病毒的影响

df_check = analysis.groupby(['date', 'country'])[['num_sold']].sum().reset_index()

• 使用 groupby 对原始数据进行分组，以日期和国家为分组依据，计算每个组的销售总量，并将结果重置索引，得到一个新的 DataFrame。

df_check['rel_gdp'] = df_check.apply(lambda s: rel_gdp_df.loc[s.country, s.date.year], axis=1)

• 为 df_check 添加一个相对 GDP 列。通过 apply 函数，根据每行的国家和年份获取相对 GDP 的值。

fig = plt.figure(figsize=(10, 4))

• 创建一个图形，设置图像大小为 (10, 4)。

for i, country in enumerate(uniques['country'][:1]):
    X = df_check[df_check['country'] == country]
    plt.plot(X['date'], X['num_sold']/X['rel_gdp'], label=country)
    plt.legend()
    plt.xlim([dt.date(2019, 10, 1), dt.date(2021, 1, 1)])

• 循环遍历国家数据中的第一个国家。
• 根据国家过滤出对应的数据子集 X。
• 使用 plt.plot 绘制归一化销售量随时间的变化，横轴为日期，纵轴为归一化销售量。
• 添加图例，显示当前国家的标签。
• 通过 plt.xlim 设置横轴范围，仅显示特定时间段。

for i in range(10, 13):
    datetime1 = dt.datetime(2019, i, 1) - dt.timedelta(days=1)
    datetime2 = dt.datetime(2019, i, 1)
    mid_datetime = datetime1 + (datetime2 - datetime1) / 2
    plt.axvline(mid_datetime, color='k', linestyle='--')

• 循环遍历从 10 到 12 的月份。
• 创建中间日期，用于在每个月份的中间绘制垂直虚线。
• 使用 plt.axvline 绘制虚线，表示每个月的分界线。

for i in range(1, 13):
    datetime1 = dt.datetime(2020, i, 1) - dt.timedelta(days=1)
    datetime2 = dt.datetime(2020, i, 1)
    mid_datetime = datetime1 + (datetime2 - datetime1) / 2
    plt.axvline(mid_datetime, color='k', linestyle='--')

• 循环遍历从 1 到 12 的月份（2020年）。
• 创建中间日期，用于在每个月份的中间绘制垂直虚线。
• 使用 plt.axvline 绘制虚线，表示每个月的分界线。
  

新冠病毒是每月效应。（这显然是一个产生的效果，因为我觉得在现实世界中，它应该是一个渐进的效果，而不是突然下降/跳跃）。将添加额外的每月功能来捕捉这种 covid 效应。

2.5分析工作日和周末的数据

绘制不同周的销售量分布。

fig = plt.figure(figsize=(10, 4))
for w in [0, 1, 2, 3]:
    w_data = analysis[analysis['week'] == w]['num_sold']
    sns.kdeplot(w_data, label=w)
plt.legend()
plt.show()

• fig = plt.figure(figsize=(10, 4)): 创建一个图形对象，设置图像大小为 (10, 4)。

• for w in [0, 1, 2, 3]:: 循环遍历周列表 [0, 1, 2, 3]。

  • w_data = analysis[analysis['week'] == w]['num_sold']: 从原始数据中选择具有指定周数的销售量数据。

  • sns.kdeplot(w_data, label=w): 使用 Seaborn 的核密度估计绘制销售量数据的分布，以及周数作为标签。

• plt.legend(): 添加图例，显示不同周数的标签。

• plt.show(): 显示第一个图形。

第二个图形的部分：

fig = plt.figure(figsize=(10, 4))
for w in [0, 1, 2, 3, 4, 5, 6]:
    w_data = analysis[analysis['week'] == w]['num_sold']
    sns.kdeplot(w_data, label=w)
plt.legend()
plt.show()

这段代码与前面的代码类似，不同之处在于循环遍历的周列表为 [0, 1, 2, 3, 4, 5, 6]，表示一周的所有七天。代码执行的步骤与前面的部分相同，绘制了第二个图形，显示不同周的销售量分布，同时添加了相应的图例。

可以看出周一至周四的分布相同，而周五至周日的分布不同

（三）线性回归（第一轮）

3.1 SMAPE指标

当面对时间序列预测问题时，有时我们需要一个用于评估预测准确性的指标。对称平均绝对百分比误差（SMAPE）是这样一种指标，它量化了预测值与实际值之间的差异，并且相对于实际值的大小进行了归一化。

下面逐步解释SMAPE的公式和含义。

SMAPE的公式是：

$\mathrm{SMAPE}=\frac{200}{n}{\sum }_{t=1}^n\frac{|F_t-A_t|}{\left(|A_t|+|F_t|\right)}\propto {\sum }_{t=1}^n\frac{|F_t/A_t-1|}{\left(1+|F_t|/|A_t|\right)}={\sum }_{F_t>A_t}1-\frac{2}{\left(1+F_t/A_t\right)}+{\sum }_{F_t<A_t}-1+\frac{2}{\left(1+F_t/A_t\right)}$

这里的符号含义为：

• $(n)$：时间序列数据的样本数量。
• $(F_t)$：预测值。
• $(A_t)$：实际值。

该公式分为两部分：分子部分是预测值与实际值的差异，分母部分是实际值和预测值的绝对值之和。这样的设计使得指标对于实际值的大小进行了标准化，避免了在评估中出现过大或过小的值。

SMAPE的取值范围在0到200之间，值越低表示预测值与实际值越接近，准确性越高。值越高则表示预测值与实际值之间的差异越大。

后面的等式意义在于，将预测值与实际值之间的差异转化为比例，并按比例进行加权。这有助于将指标解释为“偏离”的程度。

对数转换后的销售量被用作目标，因为对数变换后的值对预测偏离1的情况更为敏感。

在对数空间中，偏离1的变化更加显著，从而使模型能够更精确地捕捉预测的误差。对数变换也有助于将乘法操作转化为加法操作，使得模型更容易学习。

简单来说，这个指标是要使 𝐹𝑡/𝐴𝑡 尽可能接近1。因此，使用对数变换后的销售量作为目标，因为 𝑙𝑛$(A_t)$−𝑙𝑛$(F_t)$对于$F_t$/ $A_t$ 偏离1的敏感度比 $F_t$/-$A_t$ 更高。

3.2 模型原理:

模型描述：
该模型的目标是预测销售量（sold），其由以下几个因素相乘得到：

• 国家和年份对应的国内生产总值（GDP）
• 固定不变的商店系数（const(store)）
• 一系列包括正弦/余弦波（sine/cosine waves）的产品特征
• 假期效应（holiday）
• 周末效应（weekday）
• COVID-19疫情影响（covid）

$\text{sold}=\text{GDP}(\text{country},\text{year})\times \text{const}(\text{store})\times \left[\text{sine/cosine waves}(\text{product})+\text{holiday}+\text{weekday}+\text{covid}\right]$

取对数处理：
对上述模型的两侧取自然对数，得到以下等式：

$\ln (\text{sold})=\ln (\text{GDP})+b_1\ast \text{const(store)}+\ln \left[\text{waves(product)}+\text{holiday}+\text{weekend}+\text{covid}\right]$

对最后一项的近似：
由于波浪效应（wave）相对于其他影响较大，所以可以对最后一项进行近似，得到：

$\ln \left[\text{waves(product)}+\text{holiday}+\text{weekend}+\text{covid}\right]\approx \ln \left[\text{waves(product)}\right]+\left[\text{holiday effect}+\text{weekend}+\text{covid effect}\right]$

最终模型：
根据上述近似，最终的预测模型可以表示为：
$\ln (\text{sold})=\ln (\text{GDP})+b_1\text{const(store)}+b_2\text{(product)}+b_3\text{waves(product)}+b_4\text{holiday effect}+b_5\text{weekend}+b_6\text{covid}$

国内生产总值（GDP）作为拟合预测因子时，其系数约为1.00001，意味着其影响较小。

对于假期效应，使用多个系数的模型，以捕捉假期对销售量的影响。
每个假期会影响假期当天以及之后的一些天，因此该模型中包含一系列的系数（如 $b_{40}$、$b_{41}$ 等），分别对应于假期后的第0天、第1天等，直到假期后的第10天。

3.3 整合数据

把原始数据中的日期和销售量信息进行了加工和转换，以便于后续的分析、可视化和建模任务。把销售量取对数可能有助于将数据变换为更线性的形式，便于进行回归等建模操作。

df = pd.concat([train, test], axis=0)
df['log'] = np.log(df['num_sold'])
# add date
df['date'] = pd.to_datetime(df['date'])
df['year'] = df['date'].dt.year
df['week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df['day'] = df['date'].dt.day
df['time_no'] = (df['date'] - dt.datetime(2017, 1, 1)) // dt.timedelta(days=1)
df.loc[df['date'] > dt.datetime(2020, 2, 29), 'time_no'] -= 1

进行数据处理和特征工程，整理成适用于后续分析和建模的格式：

1. 合并数据集：将训练集和测试集按行合并，创建一个名为 df 的新数据集。
2. 添加对数列：创建一个名为 'log' 的新列，其中存储了销售量（ num_sold ）的自然对数，为了将销售量转换为更适合建模的线性尺度。
3. 处理日期信息：将原始日期（ date 列）转换为 pandas 中的日期时间格式，并从中提取出年份（ year 列）、星期几（ week 列）、月份（ month 列）和日期（ day 列）信息。
4. 计算时间编号：通过计算日期与指定日期（2017-01-01）之间的天数差，创建一个名为 'time_no' 的新列，表示日期距离起始日期经过的天数。对于2020年2月29日之后的日期，将 'time_no' 减去1，受闰年的影响。

# Define the years and countries
years = [2017, 2018, 2019, 2020, 2021, 2022]
countries = uniques['country'][:]

定义了一个包含年份（years）和国家名称（countries）的列表。年份范围从2017到2022。国家名称从之前已获取的 uniques 数据结构中获取。

# Initialize an empty list to hold DataFrames
dfs = []

创建了一个空列表 dfs，用于存储后面生成的数据框（DataFrame）。

# Generate holidays for each country and year
for year in years:
    for country in countries:
        for date, holiday_name in sorted(holidays.CountryHoliday(country, years=year).items()):
            df_0 = pd.DataFrame({"date": [date], "country": [country]})
            dfs.append(df_0)

通过嵌套的循环，对每个国家和年份生成假期信息。对于每个国家和年份，代码通过调用 holidays.CountryHoliday(country, years=year) 函数，获取了该国家在给定年份的假期信息。将每个假期的日期和国家名称构建成一个名为 df_0 的数据框，并将它添加到 dfs 列表中。

# Concatenate all the DataFrames
df_holidays = pd.concat(dfs, ignore_index=True)

使用 pd.concat 函数将所有生成的数据框从列表 dfs 连接在一起，创建了一个名为 df_holidays 的大数据框。ignore_index=True 参数表示重新生成索引。

# Convert 'date' column to datetime
df_holidays['date'] = pd.to_datetime(df_holidays['date'])

将 df_holidays 数据框中的 ‘date’ 列转换为日期时间格式，以便后续处理和分析。

df_holidays['tmp'] = 1

df_holidays 数据框添加了一个名为 ‘tmp’ 的列，并将所有行的值都设置为1。这可能是为了后续的计算或标记。

添加关于特定假期和日期的标记列，以捕捉这些特殊日期对销售量的影响

# holidays
holidays_columns = []
for i in range(0, 10):
    column = 'holiday_{}'.format(i)
    shifted = df_holidays.rename(columns={'tmp': column})
    shifted['date'] = shifted['date'] + dt.timedelta(days=i)
    df = pd.merge(df, shifted, on=['country', 'date'], how='left')
    df[column].fillna(0, inplace=True)
    df[column] = df[column]
    holidays_columns.append(column)

创建了一个空列表 holidays_columns 用于存储列名称。然后使用循环遍历0到9，为每个假期创建一个列，例如 ‘holiday_0’、‘holiday_1’ 等。
接着，通过将 df_holidays 数据框进行重命名，并将 ‘tmp’ 列重命名为上述假期列名称，得到了一个名为 shifted 的新数据框。在每次循环中，将假期的日期向后移动了 i 天，并将 ‘date’ 列进行调整。
使用 pd.merge 函数将原始数据框 df 与 shifted 数据框合并，合并键为 ‘country’ 和 ‘date’。同时，将新生成的假期列中的缺失值填充为0。
最后，将每个假期列添加到列表 holidays_columns 中。

# New Year's holiday
special_date_columns = []
for day in range(25, 32):
    column = 'day_12_{}'.format(day)
    df[column] = ((df['month'] == 12) & (df['day'] == day)).astype(float)
    special_date_columns.append(column)
for day in range(1, 11):
    column = 'day_1_{}'.format(day)
    df[column] = ((df['month'] == 1) & (df['day'] == day)).astype(float)
    special_date_columns.append(column)

首先创建了一个空列表 special_date_columns 用于存储列名称。
使用两个循环来处理特殊日期（12月25日到31日和1月1日到10日）。对于每个特殊日期，创建一个列名称，例如 ‘day_12_25’、‘day_12_26’ 等，表示12月25日、12月26日等。
通过将条件 (df['month'] == 12) & (df['day'] == day) 转换为浮点数，将特殊日期所在的行的值设为1，其余行的值设为0。最后，将每个特殊日期列添加到列表 special_date_columns 中。

获取每个国家每年的人均 GDP 数据，并进行归一化处理，生成一个关于国家和年份的 DataFrame。

def get_gdp_per_capita(country, year):
    alpha3 = {'Argentina': 'ARG', 'Canada': 'CAN',
              'Estonia': 'EST', 'Japan': 'JPN', 'Spain': 'ESP'}
    url = "https://api.worldbank.org/v2/country/{0}/indicator/NY.GDP.PCAP.CD?date={1}&format=json".format(
        alpha3[country], year)
    response = requests.get(url).json()
    return response[1][0]['value']

自定义函数 get_gdp_per_capita，用于从世界银行 API 获取指定国家和年份的人均 GDP 数据。alpha3 是一个字典，将国家名称映射到其对应的 alpha-3 代码，以便构建 API 请求的 URL。

gdp = []
for country in train.country.unique():
    row = []
    for year in range(2017, 2023):
        row.append(get_gdp_per_capita(country, year))
    gdp.append(row)

对训练数据集中的每个唯一国家循环。内部循环遍历从2017到2022年的每一年，并调用 get_gdp_per_capita 函数获取该国家每年的人均 GDP 数据。将这些数据以行的形式添加到名为 gdp 的列表中。

gdp = np.array(gdp)
gdp /= np.sum(gdp)

将 gdp 转换为 NumPy 数组，并对其进行归一化，以确保所有国家和年份的人均 GDP 数据之和为1，即转化为相对 GDP。

rel_gdp_df = pd.DataFrame(gdp, index=train.country.unique(), columns=range(2017, 2023))

将归一化后的相对 GDP 数据创建为 DataFrame，行索引为训练数据集中的唯一国家，列索引为年份。这个 DataFrame 将用于后续分析。

df['rel_gdp'] = df.apply(
    lambda s: rel_gdp_df.loc[s.country, s.date.year], axis=1)
df['log_rel_gdp'] = np.log(df['rel_gdp'])

在数据框 df 中添加关于年份的正弦波特征和产品与年份正弦波交叉特征

year_columns = ['year_sin_1', 'year_cos_1', 'year_sin_0.5', 'year_cos_0.5']
df['year_sin_1'] = np.sin(np.pi * df['time_no'] / 182.5)
df['year_cos_1'] = np.cos(np.pi * df['time_no'] / 182.5)
df['year_sin_0.5'] = np.sin(np.pi * df['time_no'] / 365.0)
df['year_cos_0.5'] = np.cos(np.pi * df['time_no'] / 365.0)

首先定义了一组年份的正弦波特征的列名。根据时间编号（即距离2017年1月1日的天数）计算了不同频率的年份正弦波和余弦波特征。这些特征将会在建模过程中用于捕捉与年份相关的周期性模式。

product_year_columns = []
for product in uniques['product']:
    df[product] = (df['product'] == product).astype(float)
    for year in year_columns:
        product_year = '{}_{}'.format(product, year)
        df[product_year] = df[product] * df[year]
        product_year_columns.append(product_year)

是为每个产品与年份正弦波交叉特征创建列。对于每个产品，将创建一个布尔列，其中元素为1表示该行的产品与当前循环的产品相同，为0表示不同。

将该布尔列与各个年份正弦波特征相乘，生成交叉特征。
这些交叉特征将在建模中用于捕捉不同产品在不同年份的周期性模式。生成的特征列名将被添加到 product_year_columns 列表中。

在数据框 df 中添加关于特定月份和年份的标志特征。

featured_month_columns = []
for month in range(3, 11):
    column = 'month_2020_{}'.format(month)
    df[column] = ((df['year'] == 2020) & (df['month'] == month)).astype(float)
    featured_month_columns.append(column)

首先定义了一组特定月份的标志特征的列名。
对于每个指定的月份（3 到 10 月），根据年份为 2020 年且月份匹配的条件，将相应的标志特征设置为 1，否则为 0。生成的标志特征列名将添加到 featured_month_columns 列表中。

each_year_columns = []
for year in [2020]:
    column = 'year_{}'.format(year)
    df[column] = ((df['year'] == year)).astype(float)
    each_year_columns.append(column)

为每个特定年份添加一个年度标志特征。对于每个指定的年份（在这里是 2020 年），将相应的年度标志特征设置为 1，表示该年份的数据行，否则为 0。生成的年度标志特征列名将添加到 each_year_columns 列表中。

在数据框 df 中添加一些关于店铺、产品和星期的标志特征。

store_columns = []
for store in uniques['store'][1:]:
    column = 'store_{}'.format(store)
    df[column] = ((df['store'] == store)).astype(float)
    store_columns.append(column)

每个店铺（除了第一个店铺）定义一个标志特征的列名。

product_columns = []
for product in uniques['product'][1:]:
    column = 'product_{}'.format(product)
    df[column] = (df['product'] == product).astype(float)
    product_columns.append(column)

每个产品（除了第一个产品）定义一个标志特征的列名。

week_columns = []
for week in range(4, 7):
    column = 'week_{}'.format(week)
    df[column] = (df['week'] == week).astype(float)
    week_columns.append(column)

每个星期（4 到 6 星期）定义一个标志特征的列名。

决定哪些列将被用作预测模型的特征

use_columns = []
use_columns.extend(special_date_columns)
use_columns.extend(product_year_columns)
use_columns.extend(holidays_columns)
use_columns.extend(week_columns)
use_columns.extend(store_columns)
use_columns.extend(product_columns)
use_columns.extend(featured_month_columns)
use_columns.extend(each_year_columns)

def fit_model(train_data, val_data, use_columns, val_year=None):
    # Separate source and target variables
    source_train = train_data[use_columns]
    target_train = train_data['log'] - \
        train_data['log_rel_gdp']  # -train_data['store_mean']

    model = linear_model.LinearRegression()
    model.fit(source_train, target_train)
    coef = pd.DataFrame(model.coef_, columns=['coef'])
    coef['column'] = source_train.columns
    coef.set_index('column', inplace=True)
    coef.loc['intercept'] = model.intercept_

    # Add predictions to the training data
    train_data['predict_log'] = model.predict(
        source_train) + train_data['log_rel_gdp']  # +train_data['store_mean']
    train_data['predict_exp'] = np.exp(train_data['predict_log'])
    train_data['smape_exp'] = 2 * (train_data['num_sold'] - train_data['predict_exp']).abs() / \
        (train_data['num_sold'] + train_data['predict_exp'])

    # Add predictions to the test data
    val_data['predict_log'] = model.predict(
        val_data[use_columns])+val_data['log_rel_gdp']  # +val_data['store_mean']
    val_data['predict_exp'] = np.exp(val_data['predict_log'])
    val_data['smape_exp'] = 2 * (val_data['num_sold'] - val_data['predict_exp']).abs() / \
        (val_data['num_sold'] + val_data['predict_exp'])

    train_mean = train_data['smape_exp'].mean()*100
    val_mean = val_data['smape_exp'].mean()*100
    val_Q1 = val_data['smape_exp'][val_data.month <= 3].mean()*100
    val_Q234 = val_data['smape_exp'][val_data.month > 3].mean()*100
    # Print the metrics
    print(
        f"train: {train_mean:.4f}, val: {val_mean:.4f}, Q1: {val_Q1:.4f}, Q234: {val_Q234:.4f}")

    return coef, model, train_mean, val_mean, val_Q1, val_Q234

定义 fit_model 的函数，用于在给定的训练数据和验证数据上执行以下步骤：

1. 将特定的列作为特征进行分离，并计算目标变量。
2. 创建一个线性回归模型，使用训练数据进行训练，然后提取模型的系数和截距。
3. 利用训练好的模型，对训练数据和验证数据进行预测，并计算预测结果的指标。
4. 计算训练数据和验证数据的平均 SMAPE（对称平均绝对百分比误差）以及在不同时间段内的 SMAPE 平均值。
5. 打印计算得到的指标，并将模型的系数、训练得分和验证得分以及其他计算结果作为输出返回。

实现了线性回归模型的训练和评估，然后将训练得分、验证得分和模型系数等信息返回给调用者。

from scipy import stats
df_used = df.copy()
# Split data into training and testing datasets
# Use the function
date_train = dt.datetime(2021, 12, 31)
train_data = df_used.loc[df_used['date'] <= date_train]
test_data = df_used.loc[(df_used['date'] >= dt.datetime(2022, 1, 1))]
from sklearn.model_selection import GroupKFold
pd.options.mode.chained_assignment = None  # default='warn'

kf = GroupKFold(n_splits=5)
train_a = []
val_a = []
val_Q1_a = []
val_Q234_a = []
for fold, (train_idx, val_idx) in enumerate(kf.split(train_data, groups=train_data.date.dt.year)):
    X_tr = train_data.iloc[train_idx]
    X_va = train_data.iloc[val_idx]
    if X_va.iloc[1].year == 2020:
        continue
    print(X_va.iloc[1].year)
    coef, model, train_mean, val_mean, val_Q1, val_Q234 = fit_model(
        X_tr, X_va, use_columns)
    train_a.append(train_mean)
    val_a.append(val_mean)
    val_Q1_a.append(val_Q1)
    val_Q234_a.append(val_Q234)

# Compute mean
train_mean = np.mean(train_a)
val_mean = np.mean(val_a)
val_Q1_mean = np.mean(val_Q1_a)
val_Q234_mean = np.mean(val_Q234_a)

# Compute standard error of the mean
train_sem = stats.sem(train_a)
val_sem = stats.sem(val_a)
val_Q1_sem = stats.sem(val_Q1_a)
val_Q234_sem = stats.sem(val_Q234_a)

print(f"Train mean: {train_mean:.4f} ± {train_sem:.4f}")
print(f"Validation mean: {val_mean:.4f} ± {val_sem:.4f}")
print(f"Validation Q1 mean: {val_Q1_mean:.4f} ± {val_Q1_sem:.4f}")
print(f"Validation Q234 mean: {val_Q234_mean:.4f} ± {val_Q234_sem:.4f}")

通过年份分为5组以进行交叉验证的方式对数据进行模型训练和验证，并输出训练和验证结果的均值及标准误差。

通过训练模型并在全部数据上进行拟合，然后应用拟合的模型对训练集数据进行预测，并计算预测结果的 SMAPE 值。

#fitting all the data
from scipy import stats
train_data = df_used.loc[df_used['date'] <= date_train]
test_data = df_used.loc[(df_used['date'] >= dt.datetime(2022, 1, 1))]
coef, model, train_mean, test_mean, test_Q1, test_Q234 = fit_model(
    train_data, test_data, use_columns)  # , alpha=150)

train_data['predict_log'] = model.predict(
    train_data[use_columns])+train_data['log_rel_gdp']
train_data['predict_exp'] = np.exp(train_data['predict_log'])
train_data['smape_exp'] = 2 * (train_data['num_sold'] - train_data['predict_exp']).abs() / \
    (train_data['num_sold'] + train_data['predict_exp'])

分析拟合模型的残差分布，通过散点图和直方图来观察残差的分布情况，并输出残差的对数标准差作为拟合模型的性能指标。

残差基本均匀。 然而，我们看到（1）：年底奇怪的直线和（2）一些奇怪的连续水平线。
可能是由于含有特殊产品；如果我们移除此产品，水平线会消失。

train_data_copy = train_data[train_data['product'] != 'Using LLMs to Win Friends and Influence People'].copy()
residuals = np.log(train_data_copy.predict_exp) - \
    np.log(train_data_copy.num_sold)
plt.figure(figsize=(18, 4))
plt.scatter(np.arange(len(residuals)), residuals, s=1)
plt.title('All residuals by row number (removing one product)')
plt.ylabel('residual')
plt.show()

绘制除以每个国家的 GDP num_sold和该特殊产品的商店

for product in ['Using LLMs to Win Friends and Influence People']:#df['product'].unique():

    countries = df['country'].unique()
    stores = df['store'].unique()
    fig, axs = plt.subplots(len(countries), len(stores), figsize=(15, 3*len(countries))) # adjust figsize as needed
    for country_idx, country in enumerate(countries):
        for store_idx, store in enumerate(stores):
            X = df[(df['country'] == country) & (df['product'] == product) & (df['store'] == store) & (df['date'] < dt.datetime(2022, 1, 1))]
            axs[country_idx, store_idx].plot(
                X['date'], X['num_sold']/X['rel_gdp'],".", label=store)
            axs[country_idx, store_idx].legend()
            if store_idx == 0:
                axs[country_idx, store_idx].set_title(country)

    plt.show()

在分析数据时，观察到Kagglazon的“num_sold”表现出良好的正弦波。但是，对于阿根廷和爱沙尼亚国家，Kaggle Learn 和 Kaggle Store 的“num_sold”显示为直线

这就提出了一个问题，即这些商店是否每天始终销售相同数量的产品。让我们深入研究数据

product='Using LLMs to Win Friends and Influence People'
country = 'Argentina'
store = 'Kaggle Learn'
fig=plt.figure(figsize=(15, 5))
X = df[(df['country'] == country) & (df['product'] == product) & (df['store'] == store) & (df['date'] < dt.datetime(2022, 1, 1))]
plt.plot(X['date'], X['num_sold'],".", label=store)
plt.title(f'{product} in {country} and {store}')
plt.show()

在残差图中，某些国家和商店的“num_sold”中观察到的水平线是“num_sold”值的每日小幅度与“num_sold”为整数的要求相结合的结果。因此，真正的基础模式（可能是正弦模式）在舍入后丢失，并且此残差无法为数据改进提供有意义的见解。

从时间序列数据中识别潜在的候选国家/地区特定假期。通过计算残差并根据国家/地区和一年中的某一天对其进行分组，它使用 z 得分检测显著偏差。这些偏差可能表明影响销售模式的假期或特殊事件

from scipy.stats import norm
train_data['dayfix'] = train_data.date.dt.dayofyear
train_data.loc[(train_data.date.dt.year != 2020) & (
    train_data.date.dt.month >= 3), 'dayfix'] += 1

print("Look for residuals beyond", norm.ppf([0.5/365, 364.5/365]))

rr = residuals.groupby([train_data.country, train_data.dayfix]).mean()
rrstd = rr.std()
print(f"Standard deviation when grouped by country and dayofyear: {rrstd:.5f}")
rrdf = pd.DataFrame({'residual': rr, 'z_score': rr / rrstd, 'date': np.datetime64(
    '2016-12-31') + pd.to_timedelta(rr.index.get_level_values(1), 'D')})
rrdf[rrdf.z_score.abs() > 3]

有很多日子有大量的残差，我检查了加拿大国家在 4 月份不同产品的预测num_sold和实际num_sold。虚线显示假日日期。可以看到，在假期的 5 天后，预测的数值远大于实际的数值。这是因为数据在此假日日没有假日效应。因此，需要从假期列表中删除此假期。其他一些假期也会发生类似的事情。
通过这种分析，我还发现日本不会在圣诞节和新年庆祝节日（2017年除外）。这就是为什么年底有大量剩余的原因。

import datetime as dt
import matplotlib.pyplot as plt

country = 'Canada'
result_graph = train_data.loc[(train_data['country'] == country) & (
    train_data['store'] == 'Kagglazon')]

products = uniques['product']
num_products = len(products)
fig, axes = plt.subplots(nrows=1, ncols=num_products, figsize=(20, 3))
i=0
year = 2017 + i
dfs = []
for date, holiday_name in sorted(holidays.CountryHoliday(country, years=year).items()):
    if date.month == 4:
        df0 = pd.DataFrame({"date": [date], "country": [country]})
        dfs.append(df0)
print(dfs)
for j, product in enumerate(products):
    for date, holiday_name in sorted(holidays.CountryHoliday(country, years=year).items()):
        if date.month == 4:
            df0 = pd.DataFrame({"date": [date], "country": [country]})
        axes[j].axvline(x=date, color='gray', linestyle='--')
    dt1, dt2 = dt.datetime(2017 + i, 4, 1), dt.datetime(2017 + i, 5, 11)
    view = result_graph.loc[(result_graph['product'] == product) & (
        result_graph['date'] > dt1)& (result_graph['date'] < dt2)]
    axes[j].plot(view['date'], view['num_sold'], color='red', label='num_sold')
    axes[j].plot(view['date'], np.around(view['predict_exp'], 0),
                 color='blue', label='predict_exp')
    axes[j].set_xlabel('time')
    axes[j].set_ylabel('num_sold')
    axes[j].set_xlim(dt1, dt2)
    axes[j].set_title(f" {product}")
    axes[j].legend()

plt.tight_layout()
plt.show()

绘制num_sold（使用 7 天平均滚动窗口视图）来分析时态模式。我们看到正弦波可以很好地拟合数据。

index = 0
view = df.loc[(df['country'] == 'Canada') &
              (df['store'] == 'Kagglazon') &
              (df['product'] == train['product'].unique()[index]),
              ['date', 'log', 'log_rel_gdp','id']]
view['log']-=view['log_rel_gdp']
for i in range(1, 7):
    shifted = view[['date', 'log']].copy()
    shifted['date'] = shifted['date'] - dt.timedelta(days=i)
    view = pd.merge(view, shifted.rename(columns={'log': 'n{}'.format(i)}),
                    how='left', on='date')
view['mean'] = view[['log', 'n1', 'n2', 'n3', 'n4', 'n5', 'n6']].mean(axis=1)
view['time_no'] = (view['date'] - dt.datetime(2017, 1, 1)
                   ) // dt.timedelta(days=1)
view.sort_values('id', inplace=True)
view = view.loc[view['date'] < dt.datetime(2019, 12, 31)]

plt.plot(view['time_no'], view['mean'])
plt.plot(view['time_no'], 9.53-0.1526*np.sin(np.pi *
         (view['time_no'] - 3) / 182.5), color='orange')
plt.axvline(x=365, color='black')
plt.axvline(x=730, color='black')
plt.title('Canada, Kagglazon, product {}'.format(
    train['product'].unique()[index]))
plt.show()

进行交叉验证，以确认某些产品只需要单个sin_1/cos_1波，某些产品只需要组合sin_0.5和cos_0.5。

(三) 线性回归（第二轮）

1. 删除一些假期
2. 新年和圣诞节不包括日本
3. 修改 SIN 和 Cos 模式

df = pd.concat([train, test], axis=0)
df['log'] = np.log(df['num_sold'])

# add date
df['date'] = pd.to_datetime(df['date'])
df['year'] = df['date'].dt.year
df['week'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df['day'] = df['date'].dt.day
df['time_no'] = (df['date'] - dt.datetime(2017, 1, 1)) // dt.timedelta(days=1)
df.loc[df['date'] > dt.datetime(2020, 2, 29), 'time_no'] -= 1

import holidays
import pandas as pd

# Define the years and countries
years = [2017, 2018, 2019, 2020, 2021, 2022]
countries = train.country.unique()

dfs = []
# Generate holidays for each country and year
for year in years:
    for country in countries:
        for date, holiday_name in sorted(holidays.CountryHoliday(country, years=year).items()):

            df_0 = pd.DataFrame({"date": [date], "country": [
                country]})
            dfs.append(df_0)

# Concatenate all the DataFrames
df_holidays = pd.concat(dfs, ignore_index=True)

# Convert 'date' column to datetime
df_holidays['date'] = pd.to_datetime(df_holidays['date'])
df_holidays['tmp'] = 1
#remove certain holidays since the data doesn't have "holiday upturn" on these holidays
df_holidays = df_holidays[~((df_holidays['date'].dt.month.isin([2,4,5,8,10])) & (df_holidays['country'] == 'Canada'))]

# it seems that Japan doesn't celebrate on this day. I remove it and it increases cross validation
df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime(
    f'2018-12-24')) & (df_holidays['country'] == 'Japan'))]


# holidays
holidays_columns = []
for i in range(0, 10):
    column = 'holiday_{}'.format(i)
    shifted = df_holidays.rename(columns={'tmp': column})
    shifted['date'] = shifted['date'] + dt.timedelta(days=i)
    df = pd.merge(df, shifted, on=['country', 'date'], how='left')
    df[column].fillna(0, inplace=True)
    df[column] = df[column]
    holidays_columns.append(column)


# New Year's holiday
special_date_columns = []
for day in range(25, 32):
    column = 'day_12_{}'.format(day)
    df[column] = ((df['month'] == 12) & (df['day'] == day) &
                  (df['country'] != 'Japan')).astype(float)
    special_date_columns.append(column)
for day in range(1, 11):
    column = 'day_1_{}'.format(day)
    df[column] = ((df['month'] == 1) & (df['day'] == day) &
                  ((df['year'] == 2017) | (df['country'] != 'Japan'))).astype(float)
    special_date_columns.append(column)
df['rel_gdp'] = df.apply(
    lambda s: rel_gdp_df.loc[s.country, s.date.year], axis=1)
df['log_rel_gdp'] = np.log(df['rel_gdp'])

# prodcut sine and cosine feature.
# sin wave
year_columns = ['year_sin_1', 'year_cos_1', 'year_sin_0.5', 'year_cos_0.5']
df['year_sin_1'] = np.sin(np.pi * df['time_no'] / 182.5)
df['year_cos_1'] = np.cos(np.pi * df['time_no'] / 182.5)
df['year_sin_0.5'] = np.sin(np.pi * df['time_no'] / 365.0)
df['year_cos_0.5'] = np.cos(np.pi * df['time_no'] / 365.0)

# I did cross validation and find that we only need simple sin and cos wave
product_year_columns = []
for product in train['product'].unique():
    df[product] = (df['product'] == product).astype(float)
    product_sin = '{}_sin'.format(product)
    product_cos = '{}_cos'.format(product)
    if product == 'Using LLMs to Train More LLMs' or product == 'Using LLMs to Win Friends and Influence People':
        df[product_sin] = df[product] * df['year_sin_0.5']
        df[product_cos] = df[product] * df['year_cos_0.5']

        product_year_columns.append(product_sin)
        product_year_columns.append(product_cos)
    elif product == 'Using LLMs to Write Better' or product == 'Using LLMs to Improve Your Coding':
        df[product_sin] = df[product] * df['year_sin_1']
        product_year_columns.append(product_sin)
    else:
        df[product_cos] = df[product] * df['year_cos_1']
        product_year_columns.append(product_cos)

# make month flag in 2020
featured_month_columns = []
for month in range(3, 11):
    column = 'month_2020_{}'.format(month)
    df[column] = ((df['year'] == 2020) & (df['month'] == month)).astype(float)
    featured_month_columns.append(column)

# I also add a yearly flag for 2020

each_year_columns = []
for year in [2020]:
    column = 'year_{}'.format(year)
    df[column] = ((df['year'] == year)).astype(float)
    each_year_columns.append(column)

store_columns = []
for store in uniques['store'][1:]:
    column = 'store_{}'.format(store)
    df[column] = ((df['store'] == store)).astype(float)
    # df[column] = ((df['store'] == store) & (df['year'] != 2020)).astype(float)
    store_columns.append(column)

# product
product_columns = []
for product in uniques['product'][1:]:
    column = 'product_{}'.format(product)
    df[column] = (df['product'] == product).astype(float)
    product_columns.append(column)

# week flag
week_columns = []
for week in range(4, 7):
    column = 'week_{}'.format(week)
    df[column] = (df['week'] == week).astype(float)
    week_columns.append(column)

use_columns = []
use_columns.extend(special_date_columns)
use_columns.extend(product_year_columns)
use_columns.extend(holidays_columns)
use_columns.extend(week_columns)
use_columns.extend(store_columns)
use_columns.extend(product_columns)
use_columns.extend(featured_month_columns)
use_columns.extend(each_year_columns)

from scipy import stats
df_used = df.copy()
# Split data into training and testing datasets
# Use the function
date_train = dt.datetime(2021, 12, 31)
train_data = df_used.loc[df_used['date'] <= date_train]
test_data = df_used.loc[(df_used['date'] >= dt.datetime(2022, 1, 1))]

kf = GroupKFold(n_splits=5)
train_a = []
val_a = []
val_Q1_a = []
val_Q234_a = []
for fold, (train_idx, val_idx) in enumerate(kf.split(train_data, groups=train_data.date.dt.year)):
    X_tr = train_data.iloc[train_idx]
    X_va = train_data.iloc[val_idx]
    if X_va.iloc[1].year == 2020:
        continue
    # if X_va.iloc[1].year == 2019:
    #     continue
    print(X_va.iloc[1].year)
    coef, model, train_mean, val_mean, val_Q1, val_Q234 = fit_model(
        X_tr, X_va, use_columns)
    train_a.append(train_mean)
    val_a.append(val_mean)
    val_Q1_a.append(val_Q1)
    val_Q234_a.append(val_Q234)

# Compute mean
train_mean = np.mean(train_a)
val_mean = np.mean(val_a)
val_Q1_mean = np.mean(val_Q1_a)
val_Q234_mean = np.mean(val_Q234_a)

# Compute standard error of the mean
train_sem = stats.sem(train_a)
val_sem = stats.sem(val_a)
val_Q1_sem = stats.sem(val_Q1_a)
val_Q234_sem = stats.sem(val_Q234_a)

print(f"Train mean: {train_mean:.4f} ± {train_sem:.4f}")
print(f"Validation mean: {val_mean:.4f} ± {val_sem:.4f}")
print(f"Validation Q1 mean: {val_Q1_mean:.4f} ± {val_Q1_sem:.4f}")
print(f"Validation Q234 mean: {val_Q234_mean:.4f} ± {val_Q234_sem:.4f}")

# fitting all the data
from scipy import stats
train_data = df_used.loc[df_used['date'] <= date_train]
test_data = df_used.loc[(df_used['date'] >= dt.datetime(2022, 1, 1))]
coef, model, train_mean, test_mean, test_Q1, test_Q234 = fit_model(
    train_data, test_data, use_columns)  # , alpha=150)

train_data['predict_log'] = model.predict(
    train_data[use_columns])+train_data['log_rel_gdp']
train_data['predict_exp'] = np.exp(train_data['predict_log'])
train_data['smape_exp'] = 2 * (train_data['num_sold'] - train_data['predict_exp']).abs() / \
    (train_data['num_sold'] + train_data['predict_exp'])

residuals = np.log(train_data.predict_exp) - np.log(train_data.num_sold)
plt.figure(figsize=(18, 4))
plt.scatter(np.arange(len(residuals)), residuals, s=1)
plt.title('All residuals by row number')
plt.ylabel('residual')
plt.show()

from scipy.stats import norm
train_data['dayfix'] = train_data.date.dt.dayofyear
train_data.loc[(train_data.date.dt.year != 2020) & (
    train_data.date.dt.month >= 3), 'dayfix'] += 1

print("Look for residuals beyond", norm.ppf([0.5/365, 364.5/365]))

rr = residuals.groupby([train_data.country, train_data.dayfix]).mean()
rrstd = rr.std()
print(f"Standard deviation when grouped by country and dayofyear: {rrstd:.5f}")
rrdf = pd.DataFrame({'residual': rr, 'z_score': rr / rrstd, 'date': np.datetime64(
    '2016-12-31') + pd.to_timedelta(rr.index.get_level_values(1), 'D')})
rrdf[rrdf.z_score.abs() > 3]

在观察数据后，我们发现只有几天表现出较大的残差。在分析假期后，我们发现加拿大和爱沙尼亚分别在 12.26 和 12.27 有假期。我们打算将这些假期分开处理，以考虑它们对数据的不同影响。此外，由于该模型已经包含了新年“大”假期的单独表示，我们的目标是通过从假期列表中删除 1.1 和 12.25 来避免多重共线性问题。这确保了对数据行为进行更准确和有意义的分析

分析假期系数

plt.plot(coef.loc['holiday_0':'holiday_9', 'coef'].values)
x = np.arange(0, 9, 0.1)
plt.xlabel('holiday 0 to holiday 9')
plt.ylabel('coefficient')
plt.show()

假日系数呈正态分布，平均值在假日后 4.5 天。做一条高斯曲线来拟合它。

from scipy.optimize import curve_fit
def gauss_function(x, a, sigma):
    return a*np.exp(-(x-4.5)**2/(sigma))

y = coef.loc['holiday_0':'holiday_9', 'coef'].values
x = np.linspace(0, len(y) - 1, len(y))
popt, pcov = curve_fit(gauss_function, x, y)
a,  sigma = popt
print(a, sigma)
plt.plot(coef.loc['holiday_0':'holiday_9', 'coef'].values)
x = np.arange(0, 9, 0.1)
plt.plot(x, a * np.exp(-(x-4.5)**2/sigma))
plt.xlabel('holiday 0 to holiday 9')
plt.ylabel('coefficient')
plt.show()

(三) 线性回归（第三轮）

1. 删除圣诞节和元旦的重复假期。
2. 分别处理加拿大和爱沙尼亚在 12.26 的假期。
3. 使用高斯曲线拟合假日系数 1-10

(三) 线性回归（第四轮）

1. 将岭回归系数绘制为正则化的函数，以找到最佳 alpha。
2. 尝试删除参数 （year_2020），并添加另一个参数 year_2019_Dec，它确实会影响 2017-2021 年期间的数据分数，但它不会影响公共排行榜上的分数。最后我提交了两个版本。没有这两个参数的版本得分更高。
3. 在这些配件之后，残差仍然具有非常轻微的月度模式，并且不是很均匀，如下所示。
4. 尝试添加 cci 数据后，并没有提高效果（呜呜呜）。

import scipy
def plot_unexplained(residuals, groups, labels=None, label_z_score=False, title=None):
    residuals_grouped = residuals.groupby(groups)
    means = residuals_grouped.mean()
    counts = residuals_grouped.count()
    z_score = np.sqrt(counts) * means / residuals.std()
    z_threshold = scipy.stats.norm.ppf(1 - 0.25 / len(means))
    m_threshold = z_threshold * residuals.std() / np.sqrt(counts.mean())
    outliers = np.abs(z_score) > z_threshold
    plt.figure(figsize=(17, 4))
    # plt.hlines([-z_threshold, +z_threshold] if label_z_score else [-m_threshold, +m_threshold], 0, len(means)-1, color='k')
    plt.bar(range(len(means)), z_score if label_z_score else means,
            color=outliers.apply(lambda b: 'r' if b else 'b'), width=0.6)
    if labels is not None:
        plt.xticks(ticks=range(len(means)), labels=labels)
    plt.ylabel('z score' if label_z_score else 'percent')
    plt.title(title)
    plt.show()


plot_unexplained(residuals, [train_data.date.dt.year, train_data.date.dt.month],
                 labels="JFMAMJJASOND" * 5,
                 title='Mean residuals of all 60 months')

(四) 模型提交

df = pd.concat([train, test], axis=0)
df['log'] = np.log(df['num_sold'])

# add date
df['date']  = pd.to_datetime(df['date'])
df['year']  = df['date'].dt.year
df['week']  = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df['day']   = df['date'].dt.day
df['time_no'] = (df['date'] - dt.datetime(2017, 1, 1)) // dt.timedelta(days=1)
df.loc[df['date'] > dt.datetime(2020, 2, 29), 'time_no'] -= 1

# Define the years and countries
years = [2017,2018,2019,2020,2021,2022]
countries = train.country.unique()

dfs = []
# Generate holidays for each country and year
for year in years:
    for country in countries:
        for date, holiday_name in sorted(holidays.CountryHoliday(country, years=year).items()):
            
            df_0 = pd.DataFrame({"date": [date], "country": [
                              country]})
            dfs.append(df_0)

# Concatenate all the DataFrames
df_holidays = pd.concat(dfs, ignore_index=True)

# Convert 'date' column to datetime
df_holidays['date'] = pd.to_datetime(df_holidays['date'])
df_holidays['tmp'] = 1

#remove certain holidays since the data doesn't have "holiday upturn" on these holidays
df_holidays = df_holidays[~((df_holidays['date'].dt.month.isin([2,4,5,8,10])) & (df_holidays['country'] == 'Canada'))]
# remove New Year and Christmas Holiday because I will handle them seperately
for country in ['Argentina', 'Canada', 'Estonia', 'Spain']:
    for year in years:
        df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime(
            f'{year}-01-01')) & (df_holidays['country'] == country))]
df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime(
    '2017-01-02')) & (df_holidays['country'] == 'Spain'))]

for country in ['Argentina', 'Canada', 'Estonia', 'Spain']:
    for year in years:
        df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime(
            f'{year}-12-25')) & (df_holidays['country'] == country))]
df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime(
    '2022-12-26')) & (df_holidays['country'] == 'Spain'))]

# Canada and Estonia has additional holiday on 26th following Christmas. I will handle them separately
for country in ['Canada', 'Estonia']:
    for year in years:
        df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime(
            f'{year}-12-26')) & (df_holidays['country'] == country))]

#Canada has additional holiday on 27, 28. I remove them and it increases cross validation
df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime('2020-12-28')) & (df_holidays['country'] == 'Canada'))]
df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime('2021-12-27')) & (df_holidays['country'] == 'Canada'))]
df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime('2021-12-28')) & (df_holidays['country'] == 'Canada'))]
df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime('2022-12-27')) & (df_holidays['country'] == 'Canada'))]

#it seems that Japan doesn't celebrate on this day. I remove it and it increases cross validation
df_holidays = df_holidays[~((df_holidays['date'] == pd.to_datetime(
            f'2018-12-24')) & (df_holidays['country'] == 'Japan'))]

def get_gdp_per_capita(country, year):
    alpha3 = {'Argentina': 'ARG', 'Canada': 'CAN',
              'Estonia': 'EST', 'Japan': 'JPN', 'Spain': 'ESP'}
    url = "https://api.worldbank.org/v2/country/{0}/indicator/NY.GDP.PCAP.CD?date={1}&format=json".format(
        alpha3[country], year)
    response = requests.get(url).json()
    return response[1][0]['value']


gdp = []
for country in train.country.unique():
    row = []
    for year in range(2017, 2023):
        row.append(get_gdp_per_capita(country, year))
    gdp.append(row)

gdp = np.array(gdp)
gdp /= np.sum(gdp)

rel_gdp_df = pd.DataFrame(gdp, index=train.country.unique(), columns=range(2017, 2023))
df['rel_gdp'] = df.apply(lambda s: rel_gdp_df.loc[s.country, s.date.year], axis=1)
df['log_rel_gdp'] = np.log(df['rel_gdp'])
# holiday
holiday_diff = [np.exp(-(i - 4.5) ** 2 / 8.5) for i in range(11)]
holidays_columns = ['holiday']
df['holiday'] = 0
for day_no, diff in enumerate(holiday_diff):
    shifted = df_holidays.copy()
    shifted['date'] = shifted['date'] + dt.timedelta(days=day_no)
    df = pd.merge(df, shifted, on=['country', 'date'], how='left')
    df['tmp'].fillna(0, inplace=True)
    df['holiday'] += df['tmp'] * diff
    df.drop(columns='tmp', inplace=True)
# New Year's holiday
special_date_columns = []
for day in range(25, 32):
    column = 'day_12_{}'.format(day)
    df[column] = ((df['month'] == 12) & (df['day'] == day) &
                  (df['country'] != 'Japan')).astype(float)
    special_date_columns.append(column)
for day in range(1, 11):
    column = 'day_1_{}'.format(day)
    df[column] = ((df['month'] == 1) & (df['day'] == day) & 
                  ((df['year'] == 2017) | (df['country'] != 'Japan'))).astype(float)
    special_date_columns.append(column)
holiday_countries = ['Estonia','Canada']
column = 'holiday_1226'
holiday_diff = [np.exp(-(i - 4.5) ** 2 / 8.5) for i in range(11)]
df[column] = 0
for day_no, diff in enumerate(holiday_diff):
    df.loc[(df['country'].isin(holiday_countries)) & (df['month'] == 12) & (df['day'] == 26 + day_no), column] = diff
    df.loc[(df['country'].isin(holiday_countries)) & (df['month'] ==  1) & (df['day'] == -5 + day_no), column] = diff

special_date_columns.append(column)

# sin wave
year_columns = ['year_sin_1', 'year_cos_1', 'year_sin_0.5', 'year_cos_0.5']
df['year_sin_1']   = np.sin(np.pi * df['time_no'] / 182.5)
df['year_cos_1']   = np.cos(np.pi * df['time_no'] / 182.5)
df['year_sin_0.5'] = np.sin(np.pi * df['time_no'] / 365.0)
df['year_cos_0.5'] = np.cos(np.pi * df['time_no'] / 365.0)

# prodcut feature.
#I did crossvalidation and find that we only need simple sin and cos wave
product_year_columns = []
for product in train['product'].unique():
    df[product] = (df['product'] == product).astype(float)
    product_sin = '{}_sin'.format(product)
    product_cos = '{}_cos'.format(product)
    if product == 'Using LLMs to Train More LLMs' or product == 'Using LLMs to Win Friends and Influence People':
        df[product_sin] = df[product] * df['year_sin_0.5']
        df[product_cos] = df[product] * df['year_cos_0.5']

        product_year_columns.append(product_sin)
        product_year_columns.append(product_cos)
    elif product == 'Using LLMs to Write Better' or product== 'Using LLMs to Improve Your Coding':
        df[product_sin] = df[product] * df['year_sin_1']
        product_year_columns.append(product_sin)    
    else:
        df[product_cos] = df[product] * df['year_cos_1']
        product_year_columns.append(product_cos)


# make month flag in 2020
featured_month_columns = []
for month in range(3, 11):
    column = 'month_2020_{}'.format(month)
    df[column] = ((df['year'] == 2020) & (df['month'] == month)).astype(float)
    featured_month_columns.append(column)
# week flag
week_columns = []
for week in range(4, 7):
    column = 'week_{}'.format(week)
    df[column] = (df['week'] == week).astype(float)
    week_columns.append(column)
store_columns = []

for store in train.store.unique()[1:]:
    column = 'store_{}'.format(store)
    df[column] = ((df['store'] == store)).astype(float)
    #df[column] = ((df['store'] == store) & (df['year'] != 2020)).astype(float)
    store_columns.append(column)
# product
product_columns = []
for product in train['product'].unique()[1:]:
    column = 'product_{}'.format(product)
    df[column] = (df['product'] == product).astype(float)
    product_columns.append(column)


from sklearn.linear_model import Ridge
# decide use columns
use_columns = []
use_columns.extend(special_date_columns)
use_columns.extend(product_year_columns)
use_columns.extend(holidays_columns)
use_columns.extend(week_columns)
use_columns.extend(store_columns)
use_columns.extend(product_columns)
use_columns.extend(featured_month_columns)


# learning
df_used = df.copy()
date = dt.datetime(2021, 12, 31)
df_used = df_used.loc[df_used['date'] <= date]
source = df_used[use_columns]
target = df_used['log']-df_used['log_rel_gdp']

model = Pipeline([
    ('standardize', StandardScaler()),
    ('linear_reg', Ridge(alpha=150, tol=0.00001, max_iter=10000))

])
model.fit(source, target)
linear_reg_model = model.named_steps['linear_reg']
coef = pd.DataFrame(linear_reg_model.coef_, columns=['coef'])
coef['column'] = source.columns
coef.set_index('column', inplace=True)
coef.loc['intercept'] = linear_reg_model.intercept_


coef

# show results
df['predict_log'] = model.predict(df[use_columns])+df['log_rel_gdp']
df['predict_exp'] = np.exp(df['predict_log'])


df['smape_log'] = 2 * (df['log'] - df['predict_log']).abs() / (df['log'] + df['predict_log'])

df['smape_exp'] = 2 * (df['num_sold'] - df['predict_exp']).abs() / (df['num_sold'] + df['predict_exp'])

result = df.loc[df['date'] <= date]
print('e smape  = {}'.format(result['smape_exp'].mean()))

(五) 模型预测

result = df.sort_values('id')
result_2 = result.copy()
result_2.loc[(result_2['year'] == 2022)&(result_2['country'] == 'Argentina'), 'predict_exp'] *= 3.372
result_2.loc[(result_2['year'] == 2022)&(result_2['country'] == 'Spain'), 'predict_exp'] *= 1.6
result_2.loc[(result_2['year'] == 2022)&(result_2['country'] == 'Japan'), 'predict_exp'] *= 1.394
result_2.loc[(result_2['year'] == 2022)&(result_2['country'] == 'Estonia'), 'predict_exp'] *= 1.651
result_2.loc[(result_2['year'] == 2022)&(result_2['country'] == 'Canada'), 'predict_exp'] *= 0.850


# First loop over products
for i, product in enumerate(train['product'].unique()):
    print(product)
    fig = plt.figure(figsize=(18, 3))

    # Then loop over countries for each product
    for country in train['country'].unique():
        result_graph = result_2.loc[(result['country'] == country) &
                                  (result['store'] == 'Kaggle Learn') & (result['date'] > dt.datetime(2021, 10, 31)) & (result['date'] < dt.datetime(2022, 5, 31))]
        view = result_graph.loc[result_graph['product'] == product]
        plt.plot(view['date'], view['predict_exp'],label=country)
    plt.legend()
    plt
    plt.show()