<jupyter_start><jupyter_text>載入資料<jupyter_code>!gdown --id '1rKgQnFdmzUEnVMHog75myo_oJ8zFSzr1' --output 'cabbage.xlsx'
!gdown --id '1-2GdInJS8I5QOO3nJCKC9Bc6YG2j6o5m' --output 'weather.csv'<jupyter_output><empty_output><jupyter_text>安裝套件<jupyter_code>%%capture
!pip install prophet
!pip install pmdarima
import pandas as pd
import numpy as np
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.statespace.sarimax import SARIMAX
import pmdarima as pm
from prophet import Prophet
%%capture
!wget -O TaipeiSansTCBeta-Regular.ttf https://drive.google.com/uc?id=1eGAsTN1HBpJAkeVM57_C7ccp7hbgSz3_&export=download
import matplotlib as mpl
import matplotlib.pyplot as plt 
from matplotlib.font_manager import fontManager

fontManager.addfont('TaipeiSansTCBeta-Regular.ttf')
mpl.rc('font', family='Taipei Sans TC Beta')<jupyter_output><empty_output><jupyter_text>檢視資料  甘藍<jupyter_code>df_c = pd.read_excel('cabbage.xlsx', index_col='date', parse_dates=True)
display(df_c.tail())
df_c.price.plot()
plt.title('每日甘藍批發價格')
plt.ylabel('元/公斤')
plt.xlabel('')
plt.show()
df_c['pq'] = df_c.price * df_c.quantity
df_c = df_c.resample('W').sum()
df_c.price = df_c.pq / df_c.quantity
df_c = df_c.drop(columns=['quantity', 'pq'])
display(df_c.tail())
df_c.plot()
plt.title('每週甘藍批發價格')
plt.ylabel('元/公斤')
plt.xlabel('')
plt.show()<jupyter_output><empty_output><jupyter_text>天氣<jupyter_code>cols = ['ObsTime', 'Temperature', 'Precp']
df_w = pd.read_csv('weather.csv', usecols=cols, index_col='ObsTime', parse_dates=True)
df_w.index.name = 'date' # ℃、mm
df_w.head()
df_w = df_w.resample('W').mean()
df_w.tail()<jupyter_output><empty_output><jupyter_text>合併<jupyter_code>df = pd.merge(df_c, df_w, how='outer', left_index=True, right_index=True)
df = df.dropna()
print(df.shape)
df.tail()
# df['price'] = np.log(df.price)
df['ln_price'] = np.log(df.price)
print(df.isna().sum())<jupyter_output>price          0
Temperature    0
Precp          0
ln_price       0
dtype: int64<jupyter_text>切割資料集<jupyter_code># Validation
val_df = df[-31:]

# Training
df = df[:-31]<jupyter_output><empty_output><jupyter_text>資料分析  Visualization<jupyter_code>import seaborn as sns
sns.heatmap(df.corr(), annot=True, cmap="YlGnBu")
plt.show()
sns.pairplot(df)
plt.show()<jupyter_output><empty_output><jupyter_text>Decompose<jupyter_code>decomp_results = seasonal_decompose(df.ln_price)
decomp_results.plot()
plt.show()
df.ln_price.plot()
plt.title('每週甘藍價格(取自然對數)')
plt.ylabel('ln_price')
plt.show()<jupyter_output><empty_output><jupyter_text>ADF TEST<jupyter_code># print(adfuller(df.price))
print(adfuller(df.ln_price))<jupyter_output>(-5.8374322113665205, 3.8471137598010745e-07, 2, 484, {'1%': -3.4439333076538854, '5%': -2.8675301216074414, '10%': -2.56996050355167}, -131.38165806041218)<jupyter_text>ACF與PACF<jupyter_code># Create figure
fig, (ax1, ax2) = plt.subplots(2,1, figsize=(8,8))

# Make ACF plot
plot_acf(df.ln_price, lags=52, zero=False, ax=ax1)

# Make PACF plot
plot_pacf(df.ln_price, lags=52, zero=False, ax=ax2, method='ywm')
plt.show()<jupyter_output><empty_output><jupyter_text>模型  SARIMAX$SARIMAX(p, d, q) \times (P, D, Q): Seasonal + ARIMA + Exogenous$$AR(p): Y_t = \beta_0 + \beta_1Y_{t-1} + \beta_2Y_{t-2} + \dots + \beta_pY_{t-p} + u_t$$MR(q): Y_t = \theta_0 + \theta_1u_{t-1} + \theta_2u_{t-2} + \dots + \theta_qu_{t-q} + u_t$$ARIMR(p, d, q): (1 - \sum_{i=1}^p \beta_i L^i)(1-L)^dY_t = (1 - \sum_{i=1}^q \theta_i L^i)u_t$, where $L$ stands for Lag operator.<jupyter_code># result = pm.auto_arima(df.ln_price, X=df[['Temperature', 'Precp']], seasonal=True, m=12, D=1)
%%capture
model = SARIMAX(df.ln_price, exog=df[['Temperature', 'Precp']],
                order=(3,0,0), seasonal_order=(2, 1, 0, 12))
result = model.fit()
print(result.summary())
# WEEK:  SARIMAX(3, 0, 0)x(2, 1, 0, 12) 
# MONTH: SARIMAX(2, 0, 0)x(2, 1, 0, 12)<jupyter_output>SARIMAX Results                                      
==========================================================================================
Dep. Variable:                           ln_price   No. Observations:                  487
Model:             SARIMAX(3, 0, 0)x(2, 1, 0, 12)   Log Likelihood                  -4.789
Date:                            Mon, 03 Jan 2022   AIC                             25.578
Time:                                    07:42:43   BIC                             58.884
Sample:                                01-01-2012   HQIC                            38.676
                                     - 04-25-2021                                         
Covariance Type:                              opg                                         
===============================================================================
                  coef    std err          z      P>|z|      [0.025      0.975]
---------------------[...]<jupyter_text>樣本內(訓練資料)預測<jupyter_code>predicted = result.get_prediction(start=-10)
mean = predicted.predicted_mean
conf1 = predicted.conf_int(alpha=0.32)  # 68%
conf2 = predicted.conf_int(alpha=0.05)  # 95%
mean = np.exp(predicted.predicted_mean)

conf1 = np.exp(conf1)
conf1.columns = ['lower price', 'upper price']

conf2 = np.exp(conf2)
conf2.columns = ['lower price', 'upper price']
fig, ax = plt.subplots()
df['2019':'2021'].price.plot(ax=ax, legend=False)
df['2021':].price.plot(style=' ', ax=ax, legend=False)
mean.plot()
plt.axvspan(mean.index.min(), mean.index.max(), color='grey', alpha=0.3)
plt.fill_between(conf1.index, conf1['lower price'], conf1['upper price'], color='xkcd:tomato red', facecolor='black')
plt.fill_between(conf2.index, conf2['lower price'], conf2['upper price'], color='xkcd:tomato red', alpha=0.6, facecolor='black')
plt.show()<jupyter_output><empty_output><jupyter_text>樣本外(測試資料)預測<jupyter_code>forecast = result.get_forecast(steps=len(val_df), exog=val_df[['Temperature', 'Precp']])
conf1 = forecast.conf_int(alpha=0.32)  # 68%
conf2 = forecast.conf_int(alpha=0.05)  # 95%
mean = np.exp(forecast.predicted_mean).append(df.price[-1:]).sort_index()
lastp = df[['price', 'price']][-1:]
lastp.columns = ['lower price', 'upper price']
lastp.index.name = 'ds'
conf1 = np.exp(conf1)
conf1.columns = ['lower price', 'upper price']
conf1 = conf1.append(lastp).sort_index()

conf2 = np.exp(conf2)
conf2.columns = ['lower price', 'upper price']
conf2 = conf2.append(lastp).sort_index()
fig, ax = plt.subplots()
df['2020':'2021'].price.plot(ax=ax, legend=False)
val_df.price.plot()
mean.plot()
plt.axvspan(mean.index.min(), mean.index.max(), color='grey', alpha=0.3)
plt.fill_between(conf1.index, conf1['lower price'], conf1['upper price'], color='xkcd:tomato red', facecolor='black')
plt.fill_between(conf2.index, conf2['lower price'], conf2['upper price'], color='xkcd:tomato red', alpha=0.6, facecolor='black')
ax.set_xlabel('')
ax.set_ylabel('元/公斤')
ax.set_ylim(0, 80)
plt.show()
result.plot_diagnostics() # pm.auto_arima(df, seasonal=True, m=7)
plt.show()<jupyter_output><empty_output><jupyter_text>fbporphet$y(t) = g(t) + s(t) + h(t) + e(t)$$g(t)$: trend models non-periodic changes.$s(t)$: seasonality presents periodic changes.$h(t)$: effects of holidays with irregular schedules.$e(t)$: covers idiosyncratic changes not accommodated by the model.  設定資料<jupyter_code>dfp = df.reset_index()
dfp.columns = ['ds', 'price', 'Temperature', 'Precp', 'ln_price']
dfp['y'] = dfp.ln_price
dfp.tail()
val_dfp = val_df.reset_index()
val_dfp.columns = ['ds', 'price', 'Temperature', 'Precp', 'ln_price']<jupyter_output><empty_output><jupyter_text>樣本內(訓練資料)預測<jupyter_code>m = Prophet(seasonality_mode='multiplicative', daily_seasonality=False, weekly_seasonality=True)
m.add_regressor('Temperature')
m.add_regressor('Precp')
m.fit(dfp[['ds', 'y', 'Temperature', 'Precp']])

predict = m.make_future_dataframe(periods=0, freq='W', include_history=True)
predict['Temperature'] = df[['Temperature']].reset_index(drop=True)
predict['Precp'] = df[['Precp']].reset_index(drop=True)

fcst0 = m.predict(predict)
fcst0 = fcst0.set_index(fcst0.ds, drop=True)<jupyter_output><empty_output><jupyter_text>樣本外(測試資料)預測  68%信賴區間<jupyter_code>m = Prophet(seasonality_mode='multiplicative', interval_width=0.68, daily_seasonality=False, weekly_seasonality=True)
m.add_regressor('Temperature')
m.add_regressor('Precp')
m.fit(dfp[['ds', 'y', 'Temperature', 'Precp']])

future = m.make_future_dataframe(periods=len(val_dfp), freq='W', include_history=False)
future['Temperature'] = val_dfp[['Temperature']]
future['Precp'] = val_dfp[['Precp']]

fcst1 = m.predict(future)
fcst1 = fcst1.set_index(fcst1.ds, drop=True)
# fig = m.plot(fcst)<jupyter_output><empty_output><jupyter_text>95%信賴區間<jupyter_code>m = Prophet(seasonality_mode='multiplicative', interval_width=0.95, daily_seasonality=False, weekly_seasonality=True)
m.add_regressor('Temperature')
m.add_regressor('Precp')
m.fit(dfp[['ds', 'y', 'Temperature', 'Precp']])

future = m.make_future_dataframe(periods=len(val_dfp), freq='W', include_history=False)
future['Temperature'] = val_dfp[['Temperature']]
future['Precp'] = val_dfp[['Precp']]

fcst2 = m.predict(future)
fcst2 = fcst2.set_index(fcst2.ds, drop=True)
# fig = m.plot(fcst)<jupyter_output><empty_output><jupyter_text>彙整<jupyter_code>mean = np.exp(fcst1.yhat).append(df.price[-1:]).sort_index()

conf1 = np.exp(fcst1[['yhat_lower', 'yhat_upper']])
conf1.columns = ['lower price', 'upper price']
conf1 = conf1.append(lastp).sort_index()

conf2 = np.exp(fcst2[['yhat_lower', 'yhat_upper']])
conf2.columns = ['lower price', 'upper price']
conf2 = conf2.append(lastp).sort_index()
fig, ax = plt.subplots()
df['2020':'2021'].price.plot(ax=ax, legend=False)
val_df.price.plot()
mean.plot()
plt.axvspan(mean.index.min(), mean.index.max(), color='grey', alpha=0.3)
plt.fill_between(conf1.index, conf1['lower price'], conf1['upper price'], color='xkcd:tomato red', facecolor='black')
plt.fill_between(conf2.index, conf2['lower price'], conf2['upper price'], color='xkcd:tomato red', alpha=0.6, facecolor='black')
ax.set_xlabel('')
ax.set_ylabel('元/公斤')
ax.set_ylim(0, 80)
plt.show()<jupyter_output><empty_output><jupyter_text>評估  準備評估資料  訓練資料(樣本內資料)<jupyter_code>dfins = pd.DataFrame({'y_real': df.price,
                      'y_sarimax': np.exp(result.get_prediction().predicted_mean),
                      'y_prophet': np.exp(fcst0.yhat)})
print(dfins.shape)
dfins.tail()<jupyter_output>(487, 3)<jupyter_text>測試資料(樣本外資料)<jupyter_code>dfoos = pd.DataFrame({'y_real': val_df.price,
                      'y_sarimax': np.exp(forecast.predicted_mean),
                      'y_prophet': np.exp(fcst1.yhat)})
print(dfoos.shape)
dfoos.head()<jupyter_output>(31, 3)<jupyter_text>[評估標準](https://scikit-learn.org/stable/modules/model_evaluation.htmlmean-absolute-error)<jupyter_code>from sklearn.metrics import mean_squared_error, mean_absolute_error, mean_absolute_percentage_error<jupyter_output><empty_output><jupyter_text>*   $\mbox{MSE}(Y, \widehat{Y}) = \frac{1}{T} \sum_{t=1}^T(Y_t - \widehat{Y}_t)^2$*   $\mbox{MAE}(Y, \widehat{Y}) = \frac{1}{T} \sum_{t=1}^T|Y_t - \widehat{Y}_t|$*   $\mbox{MAPE}(Y, \widehat{Y}) = \frac{1}{T} \sum_{t=1}^T\frac{|Y_t - \widehat{Y}_t|}{max(u_t, |Y_t|)}$  樣本內評估<jupyter_code>mse = [mean_squared_error(dfins.y_real, dfins.y_sarimax),
       mean_squared_error(dfins.y_real, dfins.y_prophet)]
mae = [mean_absolute_error(dfins.y_real, dfins.y_sarimax),
       mean_absolute_error(dfins.y_real, dfins.y_prophet)]
mape = [mean_absolute_percentage_error(dfins.y_real, dfins.y_sarimax),
        mean_absolute_percentage_error(dfins.y_real, dfins.y_prophet)]
dfv_ins = pd.DataFrame({'MSE': mse,
                        'MAE': mae,
                        'MAPE': mape}, index=['SARIMAX', 'prophet'])
dfv_ins
dfins.plot()
plt.show()<jupyter_output><empty_output><jupyter_text>樣本外評估<jupyter_code>mse = [mean_squared_error(dfoos.y_real, dfoos.y_sarimax),
       mean_squared_error(dfoos.y_real, dfoos.y_prophet)]
mae = [mean_absolute_error(dfoos.y_real, dfoos.y_sarimax),
       mean_absolute_error(dfoos.y_real, dfoos.y_prophet)]
mape = [mean_absolute_percentage_error(dfoos.y_real, dfoos.y_sarimax),
      mean_absolute_percentage_error(dfoos.y_real, dfoos.y_prophet)]
dfv_oos = pd.DataFrame({'MSE': mse,
                        'MAE': mae,
                        'MAPE': mape}, index=['SARIMAX', 'prophet'])
dfv_oos
dfoos.plot()
plt.show()
# for Question
mse = [mean_squared_error(np.log(dfoos.y_real), np.log(dfoos.y_sarimax)),
       mean_squared_error(np.log(dfoos.y_real), np.log(dfoos.y_prophet))]
mae = [mean_absolute_error(np.log(dfoos.y_real), np.log(dfoos.y_sarimax)),
       mean_absolute_error(np.log(dfoos.y_real), np.log(dfoos.y_prophet))]
mape = [mean_absolute_percentage_error(np.log(dfoos.y_real), np.log(dfoos.y_sarimax)),
      mean_absolute_percentage_error(np.log(dfoos.y_real), np.log(dfoos.y_prophet))]
dfv_oos = pd.DataFrame({'MSE': mse,
                        'MAE': mae,
                        'MAPE': mape}, index=['SARIMAX', 'prophet'])
dfv_oos<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18617469187782285, 'fraction_numerical': 0.032304983541299855, 'mean_word_length': 5.006436781609195, 'pattern_counts': {'":': 0, '<': 88, '<?xml version=': 0, '>': 89, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '516816', 'license': 'MIT', 'n_tokens_mistral': 5779, 'n_tokens_neox': 5461, 'n_words': 740, 'path': 'Others/Time_Series_MayToNov.ipynb', 'repo_name': 'p09323028/Agricultural-Prices', 'size': '13063'}