python-多元时间序列

该文对新西兰住宿人数数据进行时间序列分析,包括读取数据、可视化、单位根检验、差分处理、协整检验、模型建立和预测。通过对Hotel、BackPackers和Motel三类住宿数据的处理,发现数据不平稳,进行二阶差分使之平稳。然后使用VAR模型进行建模,进行模型稳定性和残差检验,最后进行了预测,预测了未来39天的住宿情况。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

数据来源于CensusAtSchool网站(https://new.censusatschool.org.nz/) 

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from statsmodels.sandbox.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.api import VAR
from statsmodels.stats.stattools import durbin_watson
from statsmodels.tsa.vector_ar.vecm import coint_johansen
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.tsa.stattools import adfuller as ADF

# 多元时间序列分析  VAR
# 1. 读取数据
dt = pd.read_csv("NZAccomodation.csv")
dt['Date']=pd.to_datetime(dt['DATE'])
data = dt.drop(['Date'],axis=1)
data.set_index('DATE',inplace=True)
print(data.describe())
print(data.info())

# 2. 可视化 图像可以看出数据具有周期性
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
fig, ax = plt.subplots(3, 1, sharex='col', sharey='row') # 分区
ax[0].plot(data['Hotel'],color='green',marker='o',linestyle='dashed',linewidth=1)
ax[1].plot(data['BackPackers'],color='blue',marker='*',linestyle='dashed',linewidth=1)
ax[2].plot(data['Motel'],color='red',marker='+',linestyle='dashed',linewidth=1)
ax[0].set_ylabel("Hotel")
ax[1].set_ylabel("BackPackers")
ax[2].set_ylabel("Motel")
ax[0].set_title('新西兰住宿人数图')
tick_spacing = 20
ax[0].xaxis.set_major_locator(ticker.MultipleLocator(tick_spacing))
plt.show()

# 3. 切分数据 75%训练 25%测试
trainnum = np.int64(data.shape[0]*0.75)
traindata = data.iloc[0:trainnum,:]
testdata = data.iloc[trainnum:data.shape[0],:]
print(traindata.shape) # (115, 3)
print(testdata.shape)  # (39, 3)

# 4. 单位根检验:检验序列平稳性
def Adf_diy(data):
    dftest = ADF(data,autolag='BIC')
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic','p-value','Lags Used','Number of Observations Used'])
    print("检验结果:")
    print(dfoutput)

Adf_diy(traindata.Hotel)        # p-value  0.211125 不平稳
Adf_diy(traindata.BackPackers)  # p-value  0.843871 不平稳
Adf_diy(traindata.Motel)        # p-value  0.983517 不平稳

# ACF(自相关图)、PACF(偏自相关图)
def Acf(data):
    f = plt.figure(facecolor='white')
    ax1 = f.add_subplot(211)
    plot_acf(data, lags=31, ax=ax1)
    ax2 = f.add_subplot(212)
    plot_pacf(data, lags=31, ax=ax2)
    plt.show()
Acf(traindata.Hotel)
Acf(traindata.BackPackers)
Acf(traindata.Motel)

# 5.差分并输出图像
data1 = traindata['Hotel']
data2 = traindata['BackPackers']
data3 = traindata['Motel']
tr_diff = traindata.diff().diff()
ts_diff1 = data1.diff().diff()   # 二阶差分
ts_diff2 = data2.diff().diff()
ts_diff3 = data3.diff().diff()
f = plt.figure(facecolor='white')
data1.plot(color='green', label='Hotel原始数据')
data2.plot(color='blue', label='BackPackers原始数据')
data3.plot(color='red', label='Motel原始数据')
ts_diff1.plot(color='green', linestyle='dashed',label='Hotel 2阶差分')
ts_diff2.plot(color='blue', linestyle='dashed',label='BackPackers 2阶差分')
ts_diff3.plot(color='red', linestyle='dashed',label='Motel 2阶差分')
# print(ts_diff1,ts_diff2,ts_diff3)
# plt.legend(loc='center left')
plt.title('原始数据(上)与差分后(下)')
plt.show()

# 检验一阶差分后序列是否平稳  二阶单整
tr_diff = tr_diff.dropna()
ts_diff1=ts_diff1.dropna()  # 去除空值
ts_diff2=ts_diff2.dropna()
ts_diff3=ts_diff3.dropna()
Adf_diy(ts_diff1)   # 一阶:p约等于0.0008 95%平稳
Adf_diy(ts_diff2)   #      p约等于0.1808 不平稳
Adf_diy(ts_diff3)   #      p约等于0.0791 90%下平稳


# 白噪声检验
def LB_test(timeseries):
    [[lb], [p]] = acorr_ljungbox(timeseries, lags=1)
    if p < 0.05:
        print(u"原始序列为非白噪声序列")
    else:
        print(u"原始序列为白噪声序列")

LB_test(ts_diff1)         # 差分后非白噪声序列
LB_test(ts_diff2)
LB_test(ts_diff3)

# 6. 协整(原始数据):存在则VECM(向量误差修正模型) 不存在则为伪回归
# EG两步协整检验法:仅适用于二元 1. 回归  2.残差单位根检验
# Johansen检验:特征根迹检验;最大特征值检验
# coint_johansen(endog, det_order, k_ar_diff)
#       det_order(int):-1-无确定性 0-常数项 1-线性趋势
#       k_ar_diff(非负int):模型中滞后差异的数量
def joh_output(res):
    output = pd.DataFrame([res.lr2, res.lr1],
                          index=['max_eig_stat', "trace_stat"]) # max_eig_stat最大特征值 trace_stat 迹检验值
    print(output.T, '\n')
    print("Critical values(90%, 95%, 99%) of max_eig_stat\n", res.cvm, '\n') # 最大特征值检验临界值
    print("Critical values(90%, 95%, 99%) of trace_stat\n", res.cvt, '\n') # 跟踪统计临界值

data = data.dropna()
joh_model1 = coint_johansen(data, 0, 1)  # k_ar_diff +1 = K
joh_output(joh_model1)  #(0,0)(0,1)(1,0)(1,1)(1,2)(1,3)(1,4)


def adjust(val, length= 6):
    return str(val).ljust(length)
def cointegration_test(df, alpha=0.05):
    """Perform Johanson's Cointegration Test and Report Summary"""
    out = coint_johansen(df,1,1)
    d = {'0.90':0, '0.95':1, '0.99':2}
    traces = out.lr1
    cvts = out.cvt[:, d[str(1-alpha)]]
    # Summary
    print('Name   :  Test Stat > C(95%)    =>   Signif  \n', '--'*20)
    for col, trace, cvt in zip(df.columns, traces, cvts):
        print(adjust(col), ': ', adjust(round(trace,2), 9), ">", adjust(cvt, 8), ' => ' , trace > cvt)
data = data.dropna()
cointegration_test(data)

# 平稳后变量选择
# 7. 阶数确定:LR;信息准则(AIC BIC SC)
#  VAR模型的滞后阶数越大,自由度就越小
#  VAR(p)  根据bic来推荐最优的模型
model = VAR(tr_diff)
for i in range(20):
    result = model.fit(i)
    print('Lag Order =', i)
    print('AIC : ', result.aic)
    print('BIC : ', result.bic, '\n')

# 参数估计 拟合模型
model_fitted = model.fit(14)
model_fitted.summary()

#8. 模型检验:模型稳定性检验;残差检验;各阶系数联合显著性检验
out = durbin_watson(model_fitted.resid)  # 杜宾-瓦特森检验
for col, val in zip(data.columns, out):
    print(adjust(col), ':', round(val, 2))  # 检验值越接近2,说明模型越好

# 9. 模型应用:格兰杰因果检验(平稳);脉冲响应;方差分解;预测
# 预测
lag_order = model_fitted.k_ar
# print(lag_order)
forecast_input = tr_diff.values[-lag_order:]
# print(forecast_input)
fc = model_fitted.forecast(y=forecast_input, steps=39)
df_forecast = pd.DataFrame(fc, index=data.index[-39:], columns=data.columns + '_2d')
# print(df_forecast)
# 还原
def invert_transformation(df_train, df_forecast,second_diff=False):
    df_fc = df_forecast.copy()
    columns = df_train.columns
    for col in columns:
        # 二阶差分还原
        if second_diff:
            df_fc[str(col) + '_1d'] = (df_train[col].iloc[-1] - df_train[col].iloc[-2]) + df_fc[
                str(col) + '_2d'].cumsum()
        # 一阶差分还原
        df_fc[str(col) + '_forecast'] = df_train[col].iloc[-1] + df_fc[str(col) + '_1d'].cumsum()
    return df_fc

df_results = invert_transformation(traindata, df_forecast,second_diff=True)
# print(df_results.Hotel_forecast,df_results.BackPackers_forecast,df_results.Motel_forecast)

# 存入csv文件
# df_results.to_csv("NZA预测.csv")

# 可视化
fig, axes = plt.subplots(3, 1, dpi=150, figsize=(10,10))
for i, (col,ax) in enumerate(zip(data.columns, axes.flatten())):
    df_results[col+'_forecast'].plot(legend=False, ax=ax).autoscale(axis='x',tight=True)
    testdata[col][-39:].plot(legend=False, ax=ax)
    ax.set_ylabel(col)
plt.show()

<think>嗯,用户想用Python和LSTM实现多元时间序列预测。我需要先理清楚步骤,然后一步步解释。首先,多元时间序列和单变量的区别在于,多元有多个特征,而单变量只有一个。所以数据准备的时候要注意特征的处理。 根据引用里的内容,用户提到的引用1和2都提到需要处理多列数据,比如训练时用[feature1, feature2,...,target],预测时需要提供相应的特征值。所以数据预处理是关键,可能需要标准化或者归一化,然后用滑动窗口生成时间步长的数据。 然后是模型构建部分。LSTM的结构需要设置输入层、隐藏层,可能还有Dropout防止过拟合。引用3里提到了PyTorch的LSTM模型,但用户可能更倾向于Keras,因为示例代码里有Keras的代码。不过两种框架都需要考虑输入形状,特别是时间步和特征数。 接下来是数据划分,训练集、验证集和测试集。要注意时间序列不能随机打乱,必须按时间顺序划分。引用2中的代码用了滑动窗口,每次取前30个时间步的数据,这可能需要用TimeSeriesSplit或者自定义方法。 模型训练时,损失函数通常用MSE或者MAE,优化器用Adam。回调函数如EarlyStopping可以防止过拟合。预测部分可能需要循环预测,尤其是多步预测的时候,需要逐步将预测值反馈回模型,如引用2中的for循环那样。评估指标方面,用户可能需要计算RMSE、MAE这些常见指标,可视化对比预测和真实值。另外,引用中提到预测时需要未来的特征值,这时候如果实际应用中无法获取未来特征,可能需要用其他方法,比如用滞后特征或者外部数据源。 可能用户的问题还包括如何处理缺失值,数据标准化应该在训练集上做然后应用到测试集,避免数据泄露。还有超参数调优,比如LSTM的层数、神经元数量、时间步长选择,这可能需要网格搜索或贝叶斯优化。最后,用户可能需要具体的代码示例,所以给出一个简单的Keras实现流程会比较好,包括数据预处理、模型构建、训练和预测的步骤。同时提醒注意事项,比如特征顺序、时间步的影响,以及多步预测的策略。</think>### 用Python和LSTM实现多元时间序列预测的步骤 #### 一、数据准备与预处理 1. **数据格式要求** 多元时间序列数据需包含多个特征列(如$X_1,X_2,...,X_n$)和1个目标列$Y$,格式为: $$[t, X_{1}, X_{2}, ..., X_{n}, Y]$$ 例如温度预测数据集可能包含湿度、气压等特征[^1]。 2. **标准化处理** 使用`MinMaxScaler`或`Z-Score`对数据进行归一化: ```python from sklearn.preprocessing import MinMaxScaler scaler = MinMaxScaler(feature_range=(0,1)) scaled_data = scaler.fit_transform(df) ``` 3. **时间步长构建** 定义滑动窗口长度(如$time\_steps=30$),将数据转换为3D格式$[样本数, 时间步长, 特征数]$: ```python X, y = [], [] for i in range(time_steps, len(scaled_data)): X.append(scaled_data[i-time_steps:i, :]) # 包含所有特征 y.append(scaled_data[i, target_column_index]) # 目标列 X_train, y_train = np.array(X), np.array(y) ``` #### 二、LSTM模型构建(Keras实现) ```python from keras.models import Sequential from keras.layers import LSTM, Dense, Dropout model = Sequential() model.add(LSTM(units=50, return_sequences=True, input_shape=(X_train.shape[1], X_train.shape[2]))) model.add(Dropout(0.2)) model.add(LSTM(units=30, return_sequences=False)) model.add(Dense(units=1)) model.compile(optimizer='adam', loss='mean_squared_error') ``` #### 三、模型训练与验证 1. **数据划分策略** 时间序列需按时间顺序划分,通常取前80%作为训练集: ```python train_size = int(len(X)*0.8) X_train, X_test = X[:train_size], X[train_size:] y_train, y_test = y[:train_size], y[train_size:] ``` 2. **训练配置** 使用`EarlyStopping`防止过拟合: ```python from keras.callbacks import EarlyStopping history = model.fit(X_train, y_train, epochs=100, batch_size=32, validation_split=0.1, callbacks=[EarlyStopping(monitor='val_loss', patience=5)]) ``` #### 四、预测与结果还原 1. **多步预测实现** 需要循环预测并将结果反馈给模型[^2]: ```python predictions = [] current_batch = X_test[0].reshape(1, time_steps, num_features) for i in range(len(X_test)): current_pred = model.predict(current_batch)[0] predictions.append(current_pred) # 更新批次数据,保留其他特征的真实值 current_batch = np.append(current_batch[:,1:,:], [np.append(X_test[i, -1, :-1], current_pred)], axis=1) ``` 2. **反标准化处理** 需构建完整的逆变换矩阵: ```python dummy_array = np.zeros((len(predictions), scaled_data.shape[1])) dummy_array[:, target_column_index] = predictions inv_predictions = scaler.inverse_transform(dummy_array)[:, target_column_index] ``` #### 五、关键注意事项 1. **特征-目标依赖关系** 预测时必须确保未来时间点的特征值可用。若实际场景无法获取未来特征,需采用滞后特征设计(如使用$t-1$时刻的特征预测$t$时刻目标)[^1]。 2. **时间步长选择** 通过自相关函数(ACF)分析确定最佳时间窗口,或使用网格搜索验证不同步长效果。 3. **多变量预测策略** 若需预测多个目标变量,修改输出层为`Dense(num_targets)`,并调整损失函数为多输出损失。
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值