import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 加载时间序列数据
# 假设df是一个包含日期和值的DataFrame
# df = pd.read_csv('your_data.csv', parse_dates=['date'], index_col='date')

# 平稳性检验
def test_stationarity(timeseries):
    result = adfuller(timeseries, autolag='AIC')
    print('ADF Statistic: %f' % result[0])
    print('p-value: %f' % result[1])
    print('Critical Values:')
    for key, value in result[4].items():
        print('\t%s: %.3f' % (key, value))

# 如果序列不平稳，进行差分
def difference(dataset, interval=1):
    diff = list()
    for i in range(interval, len(dataset)):
        value = dataset[i] - dataset[i - interval]
        diff.append(value)
    return diff

# 绘制ACF和PACF图
def plot_correlogram(timeseries):
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
    plot_acf(timeseries, lags=20, ax=ax1)
    plot_pacf(timeseries, lags=20, ax=ax2)
    plt.tight_layout()
    plt.show()

# 拟合ARIMA模型
def fit_arima(timeseries, order=(1, 1, 1)):
    model = ARIMA(timeseries, order=order)
    model_fit = model.fit()
    print(model_fit.summary())
    return model_fit

# 预测
def forecast(model_fit, steps=10):
    forecast_result = model_fit.forecast(steps=steps)
    conf_int = model_fit.get_forecast(steps=steps).conf_int()
    return forecast_result, conf_int

# 示例使用
# test_stationarity(df['value'])
# 如果不平稳，进行差分
# diff_series = difference(df['value'])
# test_stationarity(diff_series)
# 绘制ACF和PACF图确定p和q
# plot_correlogram(df['value'])
# 拟合ARIMA模型
# model_fit = fit_arima(df['value'], order=(1, 1, 1))
# 预测未来10个时间点
# forecast_result, conf_int = forecast(model_fit, steps=10)
                        

import numpy as np
from sklearn.metrics import mean_absolute_error

# 假设y_true是实际值，y_pred是预测值
y_true = np.array([10, 20, 30, 40, 50])
y_pred = np.array([12, 18, 33, 42, 48])

# 使用NumPy计算MAE
mae_numpy = np.mean(np.abs(y_true - y_pred))
print(f"使用NumPy计算的MAE: {mae_numpy}")

# 使用Scikit-learn计算MAE
mae_sklearn = mean_absolute_error(y_true, y_pred)
print(f"使用Scikit-learn计算的MAE: {mae_sklearn}")
                        

import numpy as np
from sklearn.metrics import mean_squared_error

# 假设y_true是实际值，y_pred是预测值
y_true = np.array([10, 20, 30, 40, 50])
y_pred = np.array([12, 18, 33, 42, 48])

# 使用NumPy计算RMSE
rmse_numpy = np.sqrt(np.mean((y_true - y_pred) ** 2))
print(f"使用NumPy计算的RMSE: {rmse_numpy}")

# 使用Scikit-learn计算RMSE
rmse_sklearn = np.sqrt(mean_squared_error(y_true, y_pred))
print(f"使用Scikit-learn计算的RMSE: {rmse_sklearn}")
                        

import numpy as np

# 假设y_true是实际值，y_pred是预测值
y_true = np.array([10, 20, 30, 40, 50])
y_pred = np.array([12, 18, 33, 42, 48])

# 使用NumPy计算MAPE
def calculate_mape(y_true, y_pred):
    # 避免除以零
    mask = y_true != 0
    return np.mean(np.abs((y_true[mask] - y_pred[mask]) / y_true[mask])) * 100

mape = calculate_mape(y_true, y_pred)
print(f"MAPE: {mape:.2f}%")
                        

import matplotlib.pyplot as plt
import seaborn as sns

# 假设y_true是实际值，y_pred是预测值
y_true = np.array([10, 20, 30, 40, 50])
y_pred = np.array([12, 18, 33, 42, 48])

# 创建时间索引
time_index = np.arange(len(y_true))

# 1. 预测值与实际值的对比图
plt.figure(figsize=(12, 6))
plt.plot(time_index, y_true, 'o-', label='实际值')
plt.plot(time_index, y_pred, 's--', label='预测值')
plt.xlabel('时间')
plt.ylabel('值')
plt.title('预测值与实际值的对比')
plt.legend()
plt.grid(True)
plt.show()

# 2. 残差图
residuals = y_true - y_pred
plt.figure(figsize=(12, 6))
plt.plot(time_index, residuals, 'o-')
plt.axhline(y=0, color='r', linestyle='-')
plt.xlabel('时间')
plt.ylabel('残差')
plt.title('残差图')
plt.grid(True)
plt.show()

# 3. 残差的直方图
plt.figure(figsize=(10, 6))
sns.histplot(residuals, kde=True)
plt.xlabel('残差')
plt.ylabel('频数')
plt.title('残差的分布')
plt.grid(True)
plt.show()

# 4. 预测值与实际值的散点图
plt.figure(figsize=(8, 8))
plt.scatter(y_true, y_pred)
plt.plot([min(y_true), max(y_true)], [min(y_true), max(y_true)], 'r--')
plt.xlabel('实际值')
plt.ylabel('预测值')
plt.title('预测值与实际值的散点图')
plt.grid(True)
plt.axis('equal')
plt.show()

# 5. 误差的分布图
plt.figure(figsize=(10, 6))
sns.boxplot(x=residuals)
plt.xlabel('误差')
plt.title('误差的分布')
plt.grid(True)
plt.show()
                        

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error
import warnings
warnings.filterwarnings('ignore')

# 生成示例数据
np.random.seed(42)
n = 200
ar_params = [0.7, -0.2]
ma_params = [0.4, -0.1]
residuals = np.random.normal(0, 1, n)

# 生成ARIMA(2,0,2)过程
y = np.zeros(n)
for i in range(2, n):
    y[i] = ar_params[0] * y[i-1] + ar_params[1] * y[i-2] + residuals[i] + ma_params[0] * residuals[i-1] + ma_params[1] * residuals[i-2]

# 创建时间索引
dates = pd.date_range(start='2020-01-01', periods=n, freq='D')
ts = pd.Series(y, index=dates)

# 绘制时间序列图
plt.figure(figsize=(12, 6))
plt.plot(ts)
plt.title('时间序列图')
plt.xlabel('日期')
plt.ylabel('值')
plt.grid(True)
plt.show()

# 平稳性检验
def test_stationarity(timeseries):
    result = adfuller(timeseries, autolag='AIC')
    print('ADF统计量: %f' % result[0])
    print('p值: %f' % result[1])
    print('临界值:')
    for key, value in result[4].items():
        print('\t%s: %.3f' % (key, value))

print("原始序列的平稳性检验:")
test_stationarity(ts)

# 如果序列不平稳，进行差分
# 在本例中，序列已经是平稳的，所以不需要差分
# 如果需要差分，可以使用以下代码：
# ts_diff = ts.diff().dropna()
# test_stationarity(ts_diff)
                        

# 绘制ACF和PACF图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(ts, lags=20, ax=ax1)
plot_pacf(ts, lags=20, ax=ax2)
plt.tight_layout()
plt.show()

# 根据ACF和PACF图，我们可以初步确定p和q的值
# 在本例中，PACF在2阶后截尾，ACF在2阶后截尾，因此可以尝试ARIMA(2,0,2)模型
                        

# 拟合ARIMA模型
model = ARIMA(ts, order=(2, 0, 2))
model_fit = model.fit()

# 输出模型摘要
print(model_fit.summary())

# 绘制模型诊断图
model_fit.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()

# 模型残差分析
residuals = model_fit.resid
plt.figure(figsize=(12, 6))
plt.plot(residuals)
plt.title('模型残差')
plt.xlabel('日期')
plt.ylabel('残差')
plt.axhline(y=0, color='r', linestyle='-')
plt.grid(True)
plt.show()

# 残差的ACF和PACF图
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(residuals, lags=20, ax=ax1)
plot_pacf(residuals, lags=20, ax=ax2)
plt.tight_layout()
plt.show()
                        

# 预测未来20个时间点
forecast_steps = 20
forecast_result = model_fit.get_forecast(steps=forecast_steps)
forecast = forecast_result.predicted_mean
conf_int = forecast_result.conf_int()

# 绘制预测结果
plt.figure(figsize=(12, 6))
plt.plot(ts, label='观测值')
plt.plot(forecast, color='red', label='预测值')
plt.fill_between(conf_int.index,
                 conf_int.iloc[:, 0],
                 conf_int.iloc[:, 1], color='pink', alpha=0.3, label='95%置信区间')
plt.title('ARIMA模型预测')
plt.xlabel('日期')
plt.ylabel('值')
plt.legend()
plt.grid(True)
plt.show()

# 计算预测误差（使用后20个观测值作为测试集）
train_size = n - forecast_steps
train, test = ts[:train_size], ts[train_size:]
model = ARIMA(train, order=(2, 0, 2))
model_fit = model.fit()
forecast = model_fit.forecast(steps=forecast_steps)
mse = mean_squared_error(test, forecast)
rmse = np.sqrt(mse)
print(f'预测均方误差(MSE): {mse:.4f}')
print(f'预测均方根误差(RMSE): {rmse:.4f}')
                        

# 尝试不同的p、d、q组合，选择AIC最小的模型
def optimize_arima(ts, p_range, d_range, q_range):
    best_aic = float('inf')
    best_order = None
    best_model = None
    
    for p in p_range:
        for d in d_range:
            for q in q_range:
                try:
                    model = ARIMA(ts, order=(p, d, q))
                    model_fit = model.fit()
                    aic = model_fit.aic
                    
                    if aic < best_aic:
                        best_aic = aic
                        best_order = (p, d, q)
                        best_model = model_fit
                except:
                    continue
    
    return best_order, best_model, best_aic

# 定义p、d、q的搜索范围
p_range = range(0, 4)
d_range = range(0, 2)
q_range = range(0, 4)

# 寻找最佳ARIMA模型
best_order, best_model, best_aic = optimize_arima(ts, p_range, d_range, q_range)
print(f'最佳ARIMA模型: ARIMA{best_order}')
print(f'最佳AIC值: {best_aic:.4f}')

# 输出最佳模型的摘要
print(best_model.summary())

# 使用最佳模型进行预测
forecast_steps = 20
forecast_result = best_model.get_forecast(steps=forecast_steps)
forecast = forecast_result.predicted_mean
conf_int = forecast_result.conf_int()

# 绘制预测结果
plt.figure(figsize=(12, 6))
plt.plot(ts, label='观测值')
plt.plot(forecast, color='red', label='预测值')
plt.fill_between(conf_int.index,
                 conf_int.iloc[:, 0],
                 conf_int.iloc[:, 1], color='pink', alpha=0.3, label='95%置信区间')
plt.title(f'ARIMA{best_order}模型预测')
plt.xlabel('日期')
plt.ylabel('值')
plt.legend()
plt.grid(True)
plt.show()