<jupyter_start><jupyter_text>**台区负荷短期预测**  目录1. [加载数据1](subtitle1)2. [加载数据2](subtitle2)<jupyter_code>import matplotlib.pylab as plt
%matplotlib inline
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 12, 4
import requests, pandas as pd, numpy as np
from pandas import DataFrame,Series
from io import StringIO
import time, json
from datetime import date,datetime,time
import statsmodels
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.seasonal import seasonal_decompose
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from pandas.tseries.offsets import Day
df = pd.read_csv("四川台区日负荷.csv")
df.序列 = df.序列.apply(lambda x: datetime.strptime(x,"%Y-%m-%d %H:%M:%S"))

df_series=df.set_index("序列").sort_index().值
df_series=df_series["2017-8-3":"2017-11-20"]# 11-24,25 and 12-11,12 do not have correct numbers. so dropped
df_diff=df_series.diff().dropna()
df_log = np.log(df_series)
df_log_dif = df_log.diff().dropna()
dfh_series = df_series.resample("H").mean()
dfh_dif = df_diff.resample("H").mean()
dfh_log = df_log.resample("H").mean()
dfh_log_dif = df_log_dif.resample("H").mean()
start_date=datetime(2017,10,5)
s_date=start_date.date().strftime("%Y-%m-%d")
e_date=(start_date+Day(7)).date().strftime("%Y-%m-%d")
t_date=(start_date+Day(8)).date().strftime("%Y-%m-%d")
t_date
ts = df_series[s_date:e_date]
ts_dif = df_diff[s_date:e_date]
ts_log = df_log[s_date:e_date]
ts_log_dif = df_log_dif[s_date:e_date]
#data_list = [ts, ts_dif, ts_log, ts_log_dif]
tsh = dfh_series[s_date:e_date]
tsh_dif = dfh_dif[s_date:e_date]
tsh_log = dfh_log[s_date:e_date]
tsh_log_dif = dfh_log_dif[s_date:e_date]

#ts_hours = df_series.resample("H").mean()
rcParams['figure.figsize'] = 12, 4
tsh.plot()<jupyter_output><empty_output><jupyter_text>稳定性检查<jupyter_code>def check_stationarity(timeseries,window):
    
    #Determing rolling statistics
    rolling_mean = timeseries.rolling(window=window,center=False).mean() 
    rolling_std = timeseries.rolling(window=window,center=False).std()

    #Plot rolling statistics:
    original = plt.plot(timeseries.index.to_pydatetime(), timeseries.values, color='blue',label='Original')
    mean = plt.plot(rolling_mean.index.to_pydatetime(), rolling_mean.values, color='red', label='Rolling Mean')
    std = plt.plot(rolling_std.index.to_pydatetime(), rolling_std.values, color='black', label = 'Rolling Std')
    plt.legend(loc='best')
    plt.title('Rolling Mean & Standard Deviation')
    plt.show(block=False)
    
    #Perform Dickey-Fuller test:
    print ('Results of Dickey-Fuller Test:')
    dickey_fuller_test = adfuller(timeseries, autolag='AIC')
    dfresults = pd.Series(dickey_fuller_test[0:4], index=['Test Statistic','p-value','#Lags Used','Number of Observations Used'])
    for key,value in dickey_fuller_test[4].items():
        dfresults['Critical Value (%s)'%key] = value
    print (dfresults)
check_stationarity(tsh_log,24)
#Decomposing-分解
#分解(decomposing) 可以用来把时序数据中的趋势和周期性数据都分离出来:
# 返回包含三个部分 trend（趋势部分） ， seasonal（季节性部分） 和residual (残留部分)
from statsmodels.tsa.seasonal import seasonal_decompose#分解(decomposing) 可以用来把时序数据中的趋势和周期性数据都分离出来

#df_decom=df_series[s_date:e_date]
decomposition = seasonal_decompose(tsh_log,freq=24)

trend = decomposition.trend#trend（趋势部分）
seasonal = decomposition.seasonal#seasonal（季节性部分）
residual = decomposition.resid#residual (残留部分)

#plt.subplot(411)
plt.plot(tsh_log, label='Original')
plt.legend(loc='best')
#plt.subplot(412)
plt.plot(trend, label='Trend')
plt.legend(loc='best')
#plt.subplot(413)
plt.plot(seasonal,label='Seasonality')
plt.legend(loc='best')
#plt.subplot(414)
plt.plot(residual, label='Residuals')
plt.legend(loc='best')
#plt.tight_layout()
check_stationarity(residual.dropna(),24)<jupyter_output><empty_output><jupyter_text>差分方法<jupyter_code>check_stationarity(tsh_log_dif,24)
lag_acf = acf(tsh_log_dif, nlags=24)#nlags表示要计算ACF的滞后数
lag_pacf = pacf(tsh_log_dif, nlags=24, method='ols')#ols  最小二乘法
plt.plot(lag_acf)
plt.axhline(y=0,linestyle='--',color='black')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_dif)),linestyle='--',color='black')
plt.axhline(y=1.96/np.sqrt(len(ts_log_dif)),linestyle='--',color='black')
plt.plot(lag_pacf)
plt.axhline(y=0,linestyle='--',color='black')
plt.axhline(y=-1.96/np.sqrt(len(ts_log_dif)),linestyle='--',color='black')
plt.axhline(y=1.96/np.sqrt(len(ts_log_dif)),linestyle='--',color='black')
model = ARIMA(tsh_log, order=(5, 1, 2))  
results = model.fit(disp=-1)
print(ts_log.diff().dropna()[:6])
print((ts_log-ts_log.shift()).dropna()[:6])
print(results.resid[:6])
print(results.fittedvalues[:6])
print(results.fittedvalues[:6]+results.resid[:6])
tsh_log.diff().dropna().plot()
results.fittedvalues.plot(color='red',alpha =0.3)
plt.title('RSS: %.4f  RSS1: %.4f' % (sum((results.fittedvalues-ts_log.diff().dropna())**2), sum(results.resid**2)))
sum((results.fittedvalues-tsh_log_dif).fillna(0)**2)/len(tsh_log_dif)
results.resid.plot(kind="kde")
print(results.summary())<jupyter_output>ARIMA Model Results                              
==============================================================================
Dep. Variable:                    D.值   No. Observations:                  191
Model:                 ARIMA(5, 1, 2)   Log Likelihood                  59.563
Method:                       css-mle   S.D. of innovations              0.175
Date:                Thu, 11 Oct 2018   AIC                           -101.126
Time:                        17:35:42   BIC                            -71.856
Sample:                    10-05-2017   HQIC                           -89.271
                         - 10-12-2017                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.0008      0.000     -1.623      0.1[...]<jupyter_text>In_sample prediction<jupyter_code>re_pred_0813=results.predict(datetime(2017,10,7,0,0,0),datetime(2017,10,7,23,0,0),typ="levels")
#print(results.forecast())
#results.predict(typ='levels')[e_date].plot()
re_pred_0813.plot()
dfh_log["2017-10-7"].plot()
plt.legend(loc="best")
pred_ts_log_dif = pd.Series(results.fittedvalues, copy=True)
(tsh_log.shift(1)+pred_ts_log_dif)['2017-10-7'].plot()
dfh_log["2017-10-7"].plot()
plt.legend(loc="best")<jupyter_output><empty_output><jupyter_text>Out of Sample prediction<jupyter_code>#re_pred_0813=results.predict(datetime(2017,8,13,0,0,0),datetime(2017,8,13,23,0,0),typ="levels")
#re_pred_0813.plot()
re_pred_tday=results.predict(dfh_log[t_date].index[0],dfh_log[t_date].index[-1],typ="levels")
re_pred_tday.plot()

plt.plot(results.forecast(24)[0])

(results.forecast(24)[0])==(re_pred_tday.values)
re_pred_tday_dif=results.predict(dfh_log[t_date].index[0],dfh_log[t_date].index[-1])
re_pred_tday_dif.plot()
re_pred_tday-re_pred_tday_dif
(re_pred_tday-re_pred_tday_dif).plot()<jupyter_output><empty_output><jupyter_text>*****Note:** result.predict() 和 result.fittedvalues 的值是一样的，都是预测的差分值. result.predict(typ='levels') 和 result.forecast()是一样的， 是预测的原始值（差分还原值）。而原值shift(1)以后和result.predict()相加即得到 result.predict(typ="levels")的值。<jupyter_code>results.predict().plot(alpha=0.5)
results.fittedvalues.plot(alpha =0.3)
(results.predict().values == results.fittedvalues.values).all()
results.predict(typ="levels").plot(alpha=0.3)
(tsh_log.shift(1).dropna()+results.predict()).plot(alpha=0.5)
((tsh_log.shift(1).dropna()+results.predict())==results.predict(typ="levels")).all()
results.forecast()[0]
(dfh_log.shift(1)[t_date]+re_pred_tday_dif).plot()
dfh_log[t_date].plot()
plt.legend()
train, test = tsh_log, dfh_log[t_date]
historical = [x for x in train]
predictions = list()

print('Printing Predicted vs Expected Values...')
print('\n')
for t in range(len(test)):
    model = ARIMA(historical, order=(2,1,1))
    model_fit = model.fit(disp=0)
    output = model_fit.forecast()
    yhat = output[0]
    predictions.append(float(yhat))
    observed = test[t]
    historical.append(observed)
    print('Predicted power load = %f, Observed power load = %f' % (np.exp(yhat), np.exp(observed)))

error = mean_squared_error(test, predictions)

print('\n')
print('Printing Mean Squared Error of Predictions...')
print('Test MSE: %.6f' % error)

predictions_series = pd.Series(predictions, index = test.index)
predictions_series.plot()
test.plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18896164639850327, 'fraction_numerical': 0.034845650140318055, 'mean_word_length': 5.048797736916549, 'pattern_counts': {'":': 1, '<': 24, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '409811', 'license': 'MIT', 'n_tokens_mistral': 3558, 'n_tokens_neox': 3367, 'n_words': 433, 'path': 'SCpowerload/Arima_model.ipynb', 'repo_name': 'biljiang/pyprojects', 'size': '8552'}