2021国赛C题

指导老师编写的代码，时间序列模型

第一部分代码 

import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import numpy as np
 
data1 = pd.read_excel('C_附件1.xlsx',sheet_name='企业的订货量（m³）')
data2 = pd.read_excel('C_附件1.xlsx',sheet_name='供应商的供货量（m³）')

for n in range(data1.shape[0]):
    plt.figure(figsize=[15,6])
    plt.plot(data1.iloc[n,2:],'r')
    plt.plot(data2.iloc[n,2:],'g')
    plt.title(data1.iloc[n,0])
    plt.show()

n=4
x = data1.iloc[n,2:]
y = data2.iloc[n,2:]
 
plt.figure(figsize=[15,6])
plt.plot(data1.iloc[n,2:],'r')
plt.plot(data2.iloc[n,2:],'g')
plt.title(data1.iloc[n,0])
plt.show()

data1_T = data1.T
data2_T = data2.T
 
 
 
def get_character(x,y):
    '获取重要性特征指标'
    y1 = y[x>0] #  去掉订单为零的数据
    x1 = x[x>0]
    # dy = y1-x1
    # 均值
    a = stats.trim_mean(y1,0.1)
    # 标准差
    b = y1.std()
    # 有效订货天数
    c = y1.shape[0]
    # 最大值
    d = y1.max()
    # 最小值
    f = y1.min()
    # 足额供货率
    e = y1[y1>=x1].shape[0] / y1.shape[0]
 
    out = [d,f,a,-b,c,e]
    return out
 
out = []
for n in range(data1_T.shape[1]):
    x = data1_T.iloc[2:,n] # 第n+1个供货商的订单量
    y = data2_T.iloc[2:,n] # 第n+1个供货商的供货量
#     print(get_character(x,y))
    out.append(get_character(x,y))
columns = ['最大值','最小值','截断均值','标准差','有效订货天数','足额供货率']
# data_c = pd.DataFrame(out,index=data1_T.index,columns=columns)
data_c = pd.DataFrame(out)
data_c.columns = columns
data_c.to_excel('out.xlsx')

data_c

n=267
print(out[n])
data1_T.iloc[2:,n].plot()
data2_T.iloc[2:,n].plot()

from GRA_BuildModel import GraModel
 
 
#建立灰色关联模型
 
 
data = pd.read_excel('out.xlsx')
 
model = GraModel(data, standard=True,p=0.1)
aa = model.result['cors']['value']
aa = pd.DataFrame(aa)
aaa = aa.mean(axis=1)
aac= aaa.sort_values(ascending=False)
aac

# 因子分析
from factor_analyzer import FactorAnalyzer
from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity
chi_square_value,p_value=calculate_bartlett_sphericity(data_c) # KMO检验
chi_square_value, p_value
 
# 因子分析
fa = FactorAnalyzer(25, rotation=None)
fa.fit(data_c)
# Check Eigenvalues
ev, v = fa.get_eigenvalues()
 
plt.scatter(range(1,data_c.shape[1]+1),ev)
plt.plot(range(1,data_c.shape[1]+1),ev)
plt.title('Scree Plot')
plt.xlabel('Factors')
plt.ylabel('Eigenvalue')
plt.grid()
plt.show()

fa = FactorAnalyzer(5, rotation="varimax")
fa.fit(data_c)
 
at = fa.get_factor_variance()
 
 
ya = fa.transform(data_c) * (at[1].T)
ya

p

fa.transform(data_c)

ccc=pd.DataFrame(fa.transform(data_c))
new= ccc.sum(axis=1).sort_values(ascending=False)

new

for i in range(50):
    n=new.index[i]
    print(n+1,out[n])
    plt.figure(figsize=[15,6])
    plt.plot(data1_T.iloc[2:,n],'r')
    plt.plot(data2_T.iloc[2:,n],'g')
    plt.title('i=%d,n=%d'%(i+1,n+1))
    plt.show()

第二部分代码

import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import itertools
import statsmodels.api as sm
import warnings
warnings.filterwarnings("ignore")
 
from statsmodels.tsa.stattools import adfuller  # 单位根检验
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox  # 白噪声检验
 
# 绘图支持中文
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
 
data1 = pd.read_excel('C_附件1.xlsx',sheet_name='企业的订货量（m³）')
data2 = pd.read_excel('C_附件1.xlsx',sheet_name='供应商的供货量（m³）')
 
t = pd.period_range('2017-2-1', periods=240, freq='w')
index = pd.PeriodIndex(t,freq='d')
 
data1_T = data1.T.iloc[2:,:]
data2_T = data2.T.iloc[2:,:]
data1_T.index=index
data2_T.index=index

def get_pdq(dy_24):
    '模型定阶'
    
    p = q =range(0,2)
    d = range(0,2)
    pdq = list(itertools.product(p,d,q))
    seasonal_pdq = [(x[0],x[1],x[2],24) for x in list(itertools.product(p,d,q))]
    out = []
    for param in pdq:
        for ps in seasonal_pdq:
            try:
                mod = sm.tsa.statespace.SARIMAX(dy_24,
                                               order=param,
                                               seasonal_order=ps,
                                               enforce_stationarity=False,
                                               enforce_invertibility=False)
                results = mod.fit()
#                 print('ARIMA{}*{} - AIC:{}'.format(param,ps,results.aic))
                out.append([param,ps,results.aic])
            except:
                continue
    return out

def forecast_y(par,dy):
    '预测'
    aic = pd.DataFrame([x[2] for x in par])
    i = aic.idxmin().loc[0]
    order = par[i][0]
    seasonal_order = par[i][1]
    mod = sm.tsa.statespace.SARIMAX(dy,
                                               order=order,
                                               seasonal_order=seasonal_order,
                                               enforce_stationarity=False,
                                               enforce_invertibility=False)
    results = mod.fit()
    pre = results.predict()
    forecast = results.forecast(24)
    return forecast,pre

# 预测24周的供货量
out = pd.DataFrame()
fore_out=pd.DataFrame()
for n in range(data2_T.shape[1]):#
    print('正在预测第 {} 个供货商,请稍候...'.format(n+1),end='    ')
    y = data2_T.iloc[:,n]
    dy_24 = y.diff(24).dropna().astype('float')
    
    # 模型定阶
    par = get_pdq(dy_24)
    
    aic = pd.DataFrame([x[2] for x in par])
    i = aic.idxmin().loc[0]
    order = par[i][0]
    seasonal_order = par[i][1]
    print(' --> 最优预测模型为：ARIMA{}*{}'.format(order,seasonal_order))
    
    # 向后预测24周
    fore,pre=forecast_y(par,dy_24)
    
    #还原数据
    t = pd.period_range('2021-9-12', periods=24, freq='w')
    index = pd.PeriodIndex(t,freq='d')
    fore.index = index
    tem = y[-24:].copy()
    tem.index = index
#     y = data2_T.iloc[:,n]
#     forecast =pd.concat([y,fore+tem])
#     plt.figure(figsize=[16,6])
#     plt.plot(forecast.dropna(),'g')
#     plt.plot(y,'r')
#     plt.show()
    te = fore.astype('int')+tem
    te[te<0]=0
    new_y = pd.concat([pre.astype('int'),te],axis=0)
    new_y[new_y<0]=0
    
    out['S'+str(n+1).zfill(3)] = te
    fore_out['S'+str(n+1).zfill(3)] = new_y
    
out.index = ['第'+ str(i+1).zfill(2) + '周' for i in range(24)]
out.T.to_excel('供货量预测结果_24周.xlsx')
 
fore_out.index = ['第'+ str(i+25).zfill(2) + '周' for i in range(240)]
fore_out.T.to_excel('供货量预测结果_全部.xlsx')

new_y = pd.concat([pre,te],axis=0)
 
plt.figure(figsize=[16,6])
plt.plot(y.astype('int'),'g')
plt.plot(new_y,'r')
plt.show()

forecast =pd.concat([y,fore+tem])
plt.figure(figsize=[16,6])
plt.plot(forecast.dropna(),'g')
plt.plot(y,'r')
plt.show()

y = data2_T.iloc[:,79]
dy_24 = y.diff(24).dropna().astype('float')
plt.figure(figsize=(15,6))
plot_acf(dy_24,lags=30)
plt.show()
plt.figure(figsize=(15,6))
# tsplot(cpi_d,lags=30)
plot_pacf(dy_24,lags=30)
plt.show()

# 白噪声检验
acorr_ljungbox(y.diff().diff(24).dropna().astype('float'),lags=1)

# 平稳性检验
 
adfuller(dy_24)