# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="Um6OuOdW4nEm"
# # 載入資料

# + id="HNkMiTtGXwm5"
# !gdown --id '1rKgQnFdmzUEnVMHog75myo_oJ8zFSzr1' --output 'cabbage.xlsx'
# !gdown --id '1-2GdInJS8I5QOO3nJCKC9Bc6YG2j6o5m' --output 'weather.csv'

# + [markdown] id="NBsbc7QU4qsl"
# # 安裝套件

# + id="aJ5wwhXH_7t3"
# %%capture
# !pip install prophet
# !pip install pmdarima

# + id="ck7soj8Unlqz"
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

# + id="X52mtZD_ZNQ8"
# %%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')

# + [markdown] id="0P_YGTCe4wjc"
# # 檢視資料

# + [markdown] id="qSrIq-mSYrEm"
# ## 甘藍

# + colab={"base_uri": "https://localhost:8080/", "height": 491} id="YQMUXqLlHgvy" outputId="c19a1ab9-889e-430d-e9f6-72e946fa91ed"
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()

# + id="nVykBc-An_rc"
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'])

# + colab={"base_uri": "https://localhost:8080/", "height": 502} id="fHrgtx7EyNAJ" outputId="e7939acc-4468-4ba8-e7e2-206e2c9ec4ff"
display(df_c.tail())
df_c.plot()
plt.title('每週甘藍批發價格')
plt.ylabel('元/公斤')
plt.xlabel('')
plt.show()

# + [markdown] id="5JbfSMSfY5L9"
# ## 天氣

# + colab={"base_uri": "https://localhost:8080/", "height": 238} id="-BfA3t4OY7Wg" outputId="9ed9b140-9911-4c2b-bd3c-9f41d9206900"
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()

# + colab={"base_uri": "https://localhost:8080/", "height": 238} id="uo-1LSFVy-G7" outputId="a0346f1b-7592-45a1-bdfd-fa6362701f7c"
df_w = df_w.resample('W').mean()
df_w.tail()

# + [markdown] id="64xAW1dPyP6W"
# ## 合併

# + colab={"base_uri": "https://localhost:8080/", "height": 255} id="-dnMJm0_ZAiT" outputId="71f82ae5-36c4-4467-b9d4-879cb3c66f86"
df = pd.merge(df_c, df_w, how='outer', left_index=True, right_index=True)
df = df.dropna()
print(df.shape)
df.tail()

# + id="JN4OjDcOPuoo"
# df['price'] = np.log(df.price)
df['ln_price'] = np.log(df.price)

# + colab={"base_uri": "https://localhost:8080/"} id="6wqkoOkwG9rG" outputId="8994e9a2-0025-4352-cf3e-2a93901e539e"
print(df.isna().sum())

# + [markdown] id="Bjcvhw3GRWS1"
# # 切割資料集

# + id="1kWatWTyRxZY"
# Validation
val_df = df[-31:]

# Training
df = df[:-31]

# + [markdown] id="CXZABd-tYwxc"
# # 資料分析

# + [markdown] id="NMRmVkccloH6"
# ## Visualization

# + colab={"base_uri": "https://localhost:8080/", "height": 269} id="zos1gbn_lpSU" outputId="9c0e764b-f62d-4d4b-9d5a-7de9c50b78bf"
import seaborn as sns
sns.heatmap(df.corr(), annot=True, cmap="YlGnBu")
plt.show()

# + colab={"base_uri": "https://localhost:8080/", "height": 726} id="EsmuMzHqmO0W" outputId="7a70fcb5-53ba-44c7-f111-22920ccb20e0"
sns.pairplot(df)
plt.show()

# + [markdown] id="A10AaU8VWSiJ"
# ## Decompose

# + colab={"base_uri": "https://localhost:8080/", "height": 297} id="npBhZ40wHt2d" outputId="92565180-44a2-41d6-e64f-89154c6dad95"
decomp_results = seasonal_decompose(df.ln_price)
decomp_results.plot()
plt.show()

# + colab={"base_uri": "https://localhost:8080/", "height": 295} id="1JRWzunSoNjx" outputId="a2b6c175-62a9-40ad-81c4-3bb8f6971cc6"
df.ln_price.plot()
plt.title('每週甘藍價格(取自然對數)')
plt.ylabel('ln_price')
plt.show()

# + [markdown] id="wS6M9otd486y"
# ## ADF TEST

# + colab={"base_uri": "https://localhost:8080/"} id="frPOnYbHoQHY" outputId="ce9ba1c4-5f47-4903-b5a6-7b829f7d58ba"
# print(adfuller(df.price))
print(adfuller(df.ln_price))

# + [markdown] id="7IT3qstn5D1k"
# ## ACF與PACF

# + colab={"base_uri": "https://localhost:8080/", "height": 499} id="K7R37LN1rKC0" outputId="75060ad6-c145-465a-83f1-7409df8fe158"
# 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()

# + [markdown] id="CAoG9LKSiwrN"
# # 模型

# + [markdown] id="30nuOwv15NRN"
# ## 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.

# + id="lqfisCiWshoE"
# result = pm.auto_arima(df.ln_price, X=df[['Temperature', 'Precp']], seasonal=True, m=12, D=1)

# + id="Zr9qtdA4l-kB"
# %%capture
model = SARIMAX(df.ln_price, exog=df[['Temperature', 'Precp']],
                order=(3,0,0), seasonal_order=(2, 1, 0, 12))
result = model.fit()

# + colab={"base_uri": "https://localhost:8080/"} id="WjF66VEF_XCn" outputId="a1ba840b-e7d6-4510-c958-05f5146c5cdb"
print(result.summary())
# WEEK:  SARIMAX(3, 0, 0)x(2, 1, 0, 12) 
# MONTH: SARIMAX(2, 0, 0)x(2, 1, 0, 12)

# + [markdown] id="gRBBJFchHxg6"
# ### 樣本內(訓練資料)預測

# + id="uFVI7paP9BQo"
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%

# + id="Wmubt6mefRg3"
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']

# + id="x6ejCx9N9A1w" colab={"base_uri": "https://localhost:8080/", "height": 291} outputId="a39ca972-f1d1-4bdc-8dc1-690262895641"
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()

# + [markdown] id="CPGje3hnH0Nf"
# ### 樣本外(測試資料)預測

# + id="7TynE1MJ0Zld"
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%

# + id="Nq4nfCeoHTVi"
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'

# + id="yTt0Czx2VYpP"
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()

# + colab={"base_uri": "https://localhost:8080/", "height": 281} id="bu-qlK010ZdR" outputId="22801df6-fb0b-4ca5-fbbc-8fb9b9d3d32e"
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()

# + colab={"base_uri": "https://localhost:8080/", "height": 295} id="JjsIM0zEvorj" outputId="ae704cbe-15f2-4291-c5ae-01909543ba86"
result.plot_diagnostics() # pm.auto_arima(df, seasonal=True, m=7)
plt.show()

# + [markdown] id="k0N-rpvSAs7Q"
# ## 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.

# + [markdown] id="Wdyw9Z4RJe5K"
# ### 設定資料

# + id="1918_z7HvyEr" colab={"base_uri": "https://localhost:8080/", "height": 206} outputId="cf0f0ce0-a6dd-4ece-d102-2e014197b4d8"
dfp = df.reset_index()
dfp.columns = ['ds', 'price', 'Temperature', 'Precp', 'ln_price']
dfp['y'] = dfp.ln_price
dfp.tail()

# + id="2Bw_wSuHc-yY"
val_dfp = val_df.reset_index()
val_dfp.columns = ['ds', 'price', 'Temperature', 'Precp', 'ln_price']

# + [markdown] id="CeZSvnKqm_Yg"
# ### 樣本內(訓練資料)預測

# + id="DcvmrZy1nIRr"
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)

# + [markdown] id="G29WocZcnNHF"
# ### 樣本外(測試資料)預測

# + [markdown] id="iSA1kKe3Jkvk"
# #### 68%信賴區間

# + id="0J5ABix5G4WD"
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)

# + [markdown] id="1ZsqHbnaJroI"
# #### 95%信賴區間

# + id="iKsrz5E0AvkF"
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)

# + [markdown] id="JRon2JFPnV47"
# #### 彙整

# + id="S9pvN_-wgDZJ"
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()

# + colab={"base_uri": "https://localhost:8080/", "height": 281} id="Wit6TtZ4gTMu" outputId="3f8c20d8-50e1-41f0-9055-2a13bd9a32fa"
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()

# + [markdown] id="182Sp-MjZ6Hg"
# # 評估

# + [markdown] id="NMIziSzqVftL"
# ## 準備評估資料

# + [markdown] id="acVB02WIn0yT"
# ### 訓練資料(樣本內資料)

# + colab={"base_uri": "https://localhost:8080/", "height": 224} id="KUZZACjVU8Tg" outputId="c846e61b-64c4-47d6-9482-f80881b3ac2c"
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()

# + [markdown] id="F5iZPiQB6Maf"
# ### 測試資料(樣本外資料)

# + colab={"base_uri": "https://localhost:8080/", "height": 224} id="IkwpoxAV1GqJ" outputId="71ac6d51-0c5b-4515-9c47-b9f645c83b31"
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()

# + [markdown] id="hn0FTOctWjf3"
# ## [評估標準](https://scikit-learn.org/stable/modules/model_evaluation.html#mean-absolute-error)
#

# + id="L76IffX2RIPn"
from sklearn.metrics import mean_squared_error, mean_absolute_error, mean_absolute_percentage_error

# + [markdown] id="H2KgbWpsWt8A"
# *   $\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|)}$
#

# + [markdown] id="pJtK1gAzYleQ"
# ## 樣本內評估

# + colab={"base_uri": "https://localhost:8080/", "height": 112} id="xnjzTFdXYoyh" outputId="ff12ddb2-4d15-4818-b1eb-71ce6becdb1c"
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

# + colab={"base_uri": "https://localhost:8080/", "height": 265} id="UFfp-CPBapGV" outputId="e459a3ab-af14-4831-9a9c-3be6c0dad756"
dfins.plot()
plt.show()

# + [markdown] id="yOBlmyv5Z6bC"
# ## 樣本外評估

# + colab={"base_uri": "https://localhost:8080/", "height": 112} id="43efEfKhZ2Pu" outputId="a757f954-32b6-4531-e3b7-77b1c4d266c6"
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

# + colab={"base_uri": "https://localhost:8080/", "height": 277} id="FWfmJw9-a-jk" outputId="90bc58f7-0fbd-4fa0-cc57-200b0b738c8f"
dfoos.plot()
plt.show()

# + colab={"base_uri": "https://localhost:8080/", "height": 112} id="UZLYZV4A6cd6" outputId="5e184a07-de8c-4a3e-9502-e9980684d9e3"
# 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

StarcoderdataJupyterScriptsDedupFiltered