import seaborn as sns
sns.set_style('darkgrid')
import pandas as pd
df = pd.read_csv('energy_demand.csv', 
                 parse_dates=['date'], 
                 index_col = 'date'
                )
df.tail()

df.info()

df.isna().sum()

df[(df.solar_exposure.isna() > 0) | (df.rainfall.isna() > 0) ]

b = {'Y': 1, 'N': 0, 1: 1, 0: 0}
df.school_day = df.school_day.map(b)
df.holiday = df.holiday.map(b)
df.head()

df.corr()

import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1)
fig.suptitle('Demand and Price time series')
axes[0] = df.demand.plot(ax=axes[0])
axes[0].set_ylabel('demand')
axes[1] = df.price.plot(ax=axes[1])
axes[1].set_ylabel('price (log scaled)')
axes[1].set_yscale('log')
plt.show()

# create a few extra columns
import numpy as np
price_monthly = df.demand.asfreq('M')
df['month'] = df.index.month
df['quarter'] = df.index.quarter
df['week'] = df.index.isocalendar().week
df['weekdayname'] = df.index.day_name()
df['logprice_negative'] = df.price_negative.apply(lambda x: np.log(abs(x)) if x<0 else x)

import numpy as np
sns.set(rc={"figure.figsize":(12, 20)})
fig, axes = plt.subplots(4, 1)
axes[0].set_title('Average Price distributions (log transformed) by quarter')
axes[0] = sns.boxplot(data=df, y=np.log(df.price), x=df.quarter, ax=axes[0])
axes[0].set_ylabel('Log Price')
axes[0].set_xlabel('quarter')
axes[1].set_title('Average Prices (log transformed) by month')
axes[1] = sns.boxplot(data=df, y=np.log(df.price), x=df.month, ax=axes[1])
axes[1].set_ylabel('Log Price')
axes[1].set_xlabel('month season')
weekdays = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"]
axes[2].set_title('Average Prices (log transformed) per weekday')
axes[2] = sns.boxplot(data=df, y=np.log(df.price), x=df.weekdayname, order=weekdays,ax=axes[2])
axes[2].set_ylabel('Prices (log tranformed)')
axes[2].set_xlabel('weekday')
axes[3].set_title('Negative Prices only (log transformed & sign changed)')
axes[3] = sns.boxplot(data=df, y=df.logprice_negative[df.logprice_negative != 0 ], x=df.weekdayname, order=weekdays)
axes[3].set_ylabel('Negative Prices (log transformed & sign changed)')
axes[3].set_xlabel('weekday')
plt.show()

df.logprice_negative[df.logprice_negative != 0 ].count()

# frac_neg_price

avg_frac_neg_price_weekday = df.frac_neg_price.groupby(df.weekdayname).mean()

sns.set(rc={"figure.figsize":(14, 8)})
fig = avg_frac_neg_price_weekday[weekdays].plot(kind='bar', x='weekdayname',
                                #order=weekdays
                               )
fig.set_ylabel('fraction of time with negative prices')
fig.set_xlabel('weekday')
fig.set_title('Fraction of time with negative prices per weekday')
#plt.yscale('log')
plt.show()

sns.set(rc={"figure.figsize":(13, 7)})
fig, axes = plt.subplots(1, 2)

fig = sns.scatterplot(data=df, y=np.log(df.price), x='demand', hue='holiday',
                 ax=axes[0]
                 )

fig = sns.scatterplot(data=df, y=np.log(df.price), x='demand', hue='school_day',
                 ax=axes[1]
                 )
plt.show()

# Import augmented dicky-fuller test function
from statsmodels.tsa.stattools import adfuller

# Run test
adfullertest_price = adfuller(df.price)
adfullertest_demand = adfuller(df.demand)

print(f'Price TS Adfuller test statistic {adfullertest_price[0]}')
print(f'Price TS Adfuller p-value {adfullertest_price[1]}')

print('##################')

print(f'Demand TS Adfuller test statistic {adfullertest_demand[0]}')
print(f'Demand TS Adfuller p-value {adfullertest_demand[1]}')

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from statsmodels.tsa.deterministic import CalendarFourier, DeterministicProcess
from xgboost import XGBRegressor

y = df.demand

# Create trend features
dp = DeterministicProcess(index=y.index, constant=True, order=2, drop=True)
X = dp.in_sample()  # I waived deliberately to use awkward out_of_sample function

# split the data
idx_train, idx_test = train_test_split(
    y.index, test_size=500, shuffle=False,
)
X_train, X_test = X.loc[idx_train, :], X.loc[idx_test, :]
y_train, y_test = y.loc[idx_train], y.loc[idx_test]

X_test

# Fit trend model
model = LinearRegression(fit_intercept=False)
model.fit(X_train, y_train)

# Make predictions
y_fit = pd.DataFrame(model.predict(X_train),index=y_train.index, columns=['demand'])
y_pred = pd.DataFrame(model.predict(X_test), index=y_test.index, columns=['demand'])

# Plot
plt.figure(figsize=(10,6))
axs = y_train.plot(color='0.25', subplots=True, sharex=True)
axs = y_test.plot(color='0.25', subplots=True, ax=axs)
axs = y_fit.plot(color='C0', subplots=True, ax=axs)
axs = y_pred.plot(color='C3', subplots=True, ax=axs)
for ax in axs: ax.legend([])
_ = plt.suptitle("Demand on downwards trend")
plt.show()

import numpy as np

def feature_engineer_X(X):
    X['dayofweek'] = X.index.dayofweek
    X['quarter'] = X.index.quarter
    X['month'] = X.index.month
    X['year'] = X.index.year
    X['dayofyear'] = X.index.dayofyear
    X['dayofmonth'] = X.index.day
    X['weekofyear'] = X.index.isocalendar().week.astype(np.int64)
    return X

X = feature_engineer_X(X)

# remove cols for XGB
tbd = ['const' , 'trend' , 'trend_squared']
for col in tbd:
    if col in X.columns:
        X.drop(tbd, axis=1, inplace=True)
        break
        
# Create splits
X_train, X_test = X.loc[idx_train, :], X.loc[idx_test, :]
y_train, y_test = y.loc[idx_train], y.loc[idx_test]

# Create residuals (the collection of detrended series) from the training set
y_resid = y_train - y_fit.demand # y_fit is DF
y_resid.tail()

# Train XGBoost on the residuals
xgb = XGBRegressor()
xgb.fit(X_train, y_resid)

# Add the predicted residuals onto the predicted trends
y_fit_boosted = xgb.predict(X_train) + y_fit.demand
y_pred_boosted = xgb.predict(X_test) + y_pred.demand

axs = y_train.plot(
    color='0.25', figsize=(15, 8), subplots=True, 
)
axs = y_test.plot(
    color='0.25', subplots=True, sharex=True, ax=axs,
)
axs = y_fit_boosted.plot(
    color='C0', subplots=True, sharex=True, ax=axs,
)
axs = y_pred_boosted.plot(
    color='C3', subplots=True, sharex=True, ax=axs,
)
plt.legend(['observed demand (training)', 'observed demand (test)', 'demand fitted (training)', 'demand forecast (test)'])
plt.ylabel('Demand in MWh')
plt.show()

axs = y_train[-100:].plot(
    color='0.25', figsize=(15, 8), subplots=True, sharex=True,
    #title=['BuildingMaterials', 'FoodAndBeverage'],
)
axs = y_test.plot(
    color='0.25', subplots=True, ax=axs,
)
axs = y_fit_boosted[-100:].plot(
    color='C0', subplots=True,  ax=axs,
)
axs = y_pred_boosted.plot(
    color='C3', subplots=True,  ax=axs,
)
plt.legend(['observed demand (training)', 'observed demand (test)', 'demand fitted (training)', 'demand forecast (test)'])
plt.ylabel('Demand in MWh')
plt.show()

lr = LinearRegression(fit_intercept=False)
lr.fit(X, y)
y_lr = lr.predict(X)

y_resid = y - y_lr
y_resid

xgb = XGBRegressor()
xgb.fit(X, y_resid)
y_xgb = xgb.predict(X)

y_demand_fitted = pd.DataFrame(
    y_lr + y_xgb,
    index=y.index,
    columns=['demand_fitted']
)

from xgboost import plot_importance

plot_importance(xgb, title='Feature importance for demand forecasting')
plt.show()

horizon_index_12M = pd.date_range('2020-10-07', periods=365, freq='D')
horizon_index_18M = pd.date_range('2020-10-07', periods=548, freq='D')

X_12M = pd.DataFrame(horizon_index_12M.dayofweek, index=horizon_index_12M, columns=['dayofweek'])
X_18M = pd.DataFrame(horizon_index_18M.dayofweek, index=horizon_index_18M, columns=['dayofweek'])

X_12M = feature_engineer_X(X_12M)
X_18M = feature_engineer_X(X_18M)

y_demand_forecast_12M = pd.DataFrame(lr.predict(X_12M) + xgb.predict(X_12M), index=horizon_index_12M, columns=['demand_forecast'])
y_demand_forecast_18M = pd.DataFrame(lr.predict(X_18M) + xgb.predict(X_18M), index=horizon_index_18M, columns=['demand_forecast'])

axs = y_demand_forecast_12M.plot(
    color='red', figsize=(15, 8), subplots=True, sharex=True, alpha=0.8

)
axs = y_demand_forecast_18M.plot(
    color='red', linestyle='--', subplots=True, sharex=True, alpha=0.6, ax=axs,

)
axs = y_demand_fitted[-200:].plot(
    color='blue', subplots=True, sharex=True, ax=axs,
)
axs = df.demand[-200:].plot(
    color='black', marker = 'o', subplots=True, sharex=True, ax=axs, alpha=0.5
)
#axs = y_pred_boosted.plot(
#    color='C3', subplots=True, sharex=True, ax=axs,
#)
#axs.legend([])
plt.legend(['demand forecast 12 months','demand forecast 18 months', 'demand fitted', 'actual demand'])
plt.title('Demand forecast for 12 & 18 months')
plt.ylabel('Demand in MWh')
plt.show()

from statsmodels.tsa.seasonal import seasonal_decompose

def plot_periodogram(ts, detrend='linear', ax=None):
    # https://www.kaggle.com/code/ryanholbrook/seasonality/tutorial
    from scipy.signal import periodogram
    fs = pd.Timedelta("1Y") / pd.Timedelta("1D")
    freqencies, spectrum = periodogram(
        ts,
        fs=fs,
        detrend=detrend,
        window="boxcar",
        scaling='spectrum',
    )
    if ax is None:
        _, ax = plt.subplots()
    ax.step(freqencies, spectrum, color="purple")
    ax.set_xscale("log")
    ax.set_xticks([1, 2, 4, 6, 12, 26, 52, 104])
    ax.set_xticklabels(
        [
            "Annual (1)",
            "Semiannual (2)",
            "Quarterly (4)",
            "Bimonthly (6)",
            "Monthly (12)",
            "Biweekly (26)",
            "Weekly (52)",
            "Semiweekly (104)",
        ],
        rotation=30,
    )
    ax.ticklabel_format(axis="y", style="sci", scilimits=(0, 0))
    ax.set_ylabel("Variance")
    ax.set_title("Price Periodogram")
    return ax

fig, axes = plt.subplots(3, 1, figsize=(10,10), gridspec_kw={'height_ratios': [1, 1,2]})
axes[1] = plot_periodogram(df.price, ax=axes[1])
ts_decomposed = seasonal_decompose(df['price'],  model='additive')
axes[0] = ts_decomposed.trend.plot( ax=axes[0])
axes[0].set_title("Price trend component")
axes[0].set_yscale("log")
axes[2] = sns.boxplot(data=df, y='price', showfliers=False)
axes[2].set_title("Price distribution")
plt.show()

from warnings import simplefilter

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from xgboost import XGBRegressor

simplefilter("ignore")

def make_lags(ts, lags, lead_time=1):
    return pd.concat(
        {
            f'y_lag_{i}': ts.shift(i)
            for i in range(lead_time, lags + lead_time)
        },
        axis=1)


# Four weeks of lag features
y = np.log(df.price.copy())
y_orig = y.copy()
lags = int(1)
X = make_lags(y, lags=lags).fillna(0.0)

# remove cols for XGB
tbd = ['dayofweek',	'quarter',	'month',	'dayofyear',	'dayofmonth',	'weekofyear	']
for col in tbd:
    if col in X.columns:
        X.drop(tbd, axis=1, inplace=True)
        break
        

def make_multistep_target(ts, steps):
    return pd.concat(
        {f'y_step_{i + 1}': ts.shift(-i)
         for i in range(steps)},
        axis=1)


# Eight-week forecast
y = make_multistep_target(y, steps=1).dropna()

# Shifting has created indexes that don't match. Only keep times for
# which we have both targets and features.
y, X = y.align(X, join='inner', axis=0)

## Feature Engineering

X = feature_engineer_X(X)
#X.drop("year", axis=1, inplace=True)
X['demand_modeled'] = y_demand_fitted
#X.drop('y_lag_1', axis=1,  inplace=True)

X.info()

y.head(3)

X.head(4)

# Create splits
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, shuffle=False)

#from sklearn.multioutput import MultiOutputRegressor

pp = {'colsample_bytree': 0.6, 'gamma': 0.075, 'learning_rate': 0.04, 'max_depth': 3, 'min_child_weight': 2, 'n_estimators': 200, 'subsample': 0.7}
price_model = XGBRegressor(**pp)
#price_model2 = MultiOutputRegressor(XGBRegressor())
price_model.fit(X_train, y_train)

y_fit = pd.DataFrame(price_model.predict(X_train), index=X_train.index, columns=y.columns)
y_pred = pd.DataFrame(price_model.predict(X_test), index=X_test.index, columns=y.columns)

train_rmse = mean_squared_error(y_train, y_fit, squared=False)
test_rmse = mean_squared_error(y_test, y_pred, squared=False)
print((f"Tuned model: Train RMSE: {train_rmse:.3f}\n" f"Test RMSE: {test_rmse:.3f}"))

# Set Matplotlib defaults
plt.style.use("seaborn-whitegrid")
plt.rc("figure", autolayout=True, figsize=(11, 4))
plt.rc(
    "axes",
    labelweight="bold",
    labelsize="large",
    titleweight="bold",
    titlesize=16,
    titlepad=10,
)
plot_params = dict(
    color="0.75",
    style=".-",
    markeredgecolor="0.25",
    markerfacecolor="0.25",
)
%config InlineBackend.figure_format = 'retina'


def plot_multistep(y, every=1, ax=None, palette_kwargs=None):
    palette_kwargs_ = dict(palette='husl', n_colors=16, desat=None)
    if palette_kwargs is not None:
        palette_kwargs_.update(palette_kwargs)
    palette = sns.color_palette(**palette_kwargs_)
    if ax is None:
        fig, ax = plt.subplots()
    ax.set_prop_cycle(plt.cycler('color', palette))
    #ax.set_yscale('log')
    for date, preds in y[::every].iterrows():
        preds.index = pd.period_range(start=date, periods=len(preds))
        preds.plot(ax=ax)
    return ax

train_rmse = mean_squared_error(y_train, y_fit, squared=False)
test_rmse = mean_squared_error(y_test, y_pred, squared=False)
print((f"Train RMSE: {train_rmse:.3f}\n" f"Test RMSE: {test_rmse:.3f}"))
MAE = np.mean(np.absolute(y_test.y_step_1 - y_pred.y_step_1))
MAE_Per = 100* MAE/np.mean(y_test.y_step_1)
print(MAE)
print(MAE_Per)

axs = y_train.plot(
    color='0.25', figsize=(13, 8), marker = 'o',
    sharex=True,
)
axs = y_test.plot(
    color='0.5', #subplots=True, 
    sharex=True, ax=axs, marker = 'o',
)
axs = y_fit.plot(
    color='C0', #subplots=True, 
    sharex=True, ax=axs,
)
axs = y_pred.plot(
    color='C3', #subplots=True, 
    sharex=True, ax=axs, label="modeled price"
)
axs.legend(["observed prices (train data)", "observed prices (test data)", "modeled price (train data)", "modeled price (test data)"])
plt.ylabel('prices (log transformed)')
plt.title('Observed and modeled prices')
plt.show()

axs = y_test[250:600].plot(
    color='0.25', figsize=(13, 9),
    marker = 'o',
    sharex=True,
)

axs = y_pred[250:600].plot(
    color='C3',
    linewidth=2, #linestyle = '--',
    sharex=True, ax=axs,
    label = 'Forecast'
)

"""axs = y_pred2[250:600].plot(
    color='green',
    linewidth=2, #linestyle = '--',
    sharex=True, ax=axs,
    label = 'Forecast'
)"""

axs.legend(['Actual', 'Forecast'])
plt.ylabel('prices (log transformed)')
plt.title('Zoomed in detail view of the indiviual stepwise forecast prices (test data)')
plt.show()

axs = y_test[350:600].plot(
    color='0.25', figsize=(13, 9),
    marker = 'o',
    sharex=True,
)

axs = y_pred[350:600].plot(
    color='C3',
    linewidth=2, #linestyle = '--',
    sharex=True, ax=axs
)

"""axs = y_pred2[350:600].plot(
    color='green',
    linewidth=2, linestyle = '--',
    sharex=True, ax=axs
)
"""
axs.legend(['Actual', 'Forecast', 'Forecast (unoptimized model)'])
plt.ylabel('prices (log transformed)')
plt.title('Zoomed in detail view of the indiviual stepwise forecast prices (test data)')
plt.show()

from xgboost import plot_importance
plot_importance(price_model, title='Feature importance in our price model')
plt.show()

X_12M = pd.DataFrame(horizon_index_12M.dayofweek, index=horizon_index_12M, columns=['dayofweek'])
X_18M = pd.DataFrame(horizon_index_18M.dayofweek, index=horizon_index_18M, columns=['dayofweek'])

X_12M = make_lags(X_12M, lags=lags).fillna(0.0)
X_18M = make_lags(X_18M, lags=lags).fillna(0.0)

X_12M = feature_engineer_X(X_12M)
#X_12M.drop("year", axis=1, inplace=True)
X_12M['demand'] = y_demand_forecast_12M
X_18M = feature_engineer_X(X_18M)
#X_18M.drop("year", axis=1, inplace=True)
X_18M['demand'] = y_demand_forecast_18M
#X_12M.drop('dayofweek', axis=1, inplace=True)
#X_18M.drop('dayofweek', axis=1, inplace=True)

#y_pred = pd.DataFrame(model.predict(X_12M), index=horizon_index_12M, columns=['price_forecast'])
pricetrend_adj = -0.5

y_price_forecast_12M = pd.DataFrame(price_model.predict(X_12M)+pricetrend_adj, index=horizon_index_12M, columns=['price_forecast'])
y_price_forecast_18M = pd.DataFrame(price_model.predict(X_18M)+pricetrend_adj, index=horizon_index_18M, columns=['price_forecast'])

axs = y_price_forecast_12M.plot(
    color='red', figsize=(15, 8), subplots=True, sharex=True, alpha=0.8

)
axs = y_price_forecast_18M.plot(
    color='red', linestyle='--', subplots=True, sharex=True, alpha=0.6, ax=axs,

)
#axs = y_pred[-700:].plot(
#    color='blue', subplots=True, sharex=True, ax=axs,
#)
axs = np.log(df.price[-700:]).plot(
    color='black', marker = 'o', subplots=True, sharex=True, ax=axs, alpha=0.5
)
#axs = y_pred_boosted.plot(
#    color='C3', subplots=True, sharex=True, ax=axs,
#)
#axs.legend([])
plt.legend(['Price forecast 12 months','Price forecast 18 months', 
            #'price fitted', 
            'actual price'])
plt.title('Price forecast for 12 & 18 months')
plt.ylabel('prices (log transformed)')
plt.show()

fig, axes = plt.subplots(1, 10, sharey=True, figsize=(14,8))
fig.suptitle('Price distributions comparison between actual prices (blue) and modeled prices (orange) ex outliers')
pal = sns.color_palette("PuBu", 12); pal2=sns.color_palette("OrRd", 10)
sns.boxplot(data=df['2015'],  y='price', showfliers = False, ax=axes[0], color=pal[4])
axes[0].set_title('year 2015')
sns.boxplot(data=df['2016'],  y='price', showfliers = False, ax=axes[1], color=pal[5])
axes[1].set_title('year 2016')
sns.boxplot(data=df['2017'],  y='price', showfliers = False, ax=axes[2], color=pal[6])
axes[2].set_title('year 2017')
sns.boxplot(data=df['2018'],  y='price', showfliers = False, ax=axes[3], color=pal[7])
axes[3].set_title('year 2018')
sns.boxplot(data=df['2019'],  y='price', showfliers = False, ax=axes[4], color=pal[8])
axes[4].set_title('year 2019')
sns.boxplot(data=df['2020'],  y='price', showfliers = False, ax=axes[5], color=pal[9])
axes[5].set_title('year 2020 YTD')
sns.boxplot(data=df,  y='price', showfliers = False, ax=axes[6], color=pal[10])
axes[6].set_title("years 2015-2020")
sns.boxplot(y=np.exp(y_fit.y_step_1), showfliers = False, ax=axes[7], color=pal2[3])
axes[7].set_title("modeled prices\n training")
sns.boxplot(y=np.exp(y_price_forecast_12M['price_forecast']), showfliers = False, ax=axes[8], color=pal2[5])
axes[8].set_title("forecast prices\n 12 months")
sns.boxplot(y=np.exp(y_price_forecast_18M['price_forecast']), showfliers = False,  ax=axes[9], color=pal2[5])
axes[9].set_title("forecast prices\n 18 months")
plt.show()

import random
# set up required cols
df['stored_energy'] = 0
df['revenue_negprices'] = 0
df['revenue_regular'] = 0
df['year'] = df.index.year

# function to determine daily if and how much energy is favourable to be stored and how much revenue could have been made (ex post)
def energy_storage(df, index, row, MWh):
        # sell stored energy as soon as price is somehow favourable at avg spot price
        loc = df.index.get_loc(index) # integer index
        if df.iloc[loc-1]['stored_energy'] > 0:
            df.loc[index, 'revenue_regular'] = df.iloc[loc-1]['stored_energy'] * row.price_positive
            df.loc[index, 'stored_energy'] = 0
        
        # unfavourable market prices for selling electricity
        if row.price_negative < 0:
            # A. we store our own produced PV energy free of cost or 
            # B. we sell our storage capacity for market price and earn revenue by 'buying' electricity at negative prices
            # the chance is simplified and assumed to 50%/50% suntime/nighttime
            suntime = random.choice([0, 1])
            if row.demand_neg_price > MWh * 2 & suntime == 1: # suntime situation is large and long enough to fill 70MWh storage
                df.loc[index, 'stored_energy'] = MWh
                df.loc[index, 'revenue_negprices'] = 0 # not yet sold, but we may have avoided buying spot energy storage
            elif row.demand_neg_price > MWh * 2 & suntime == 0: # night
                df.loc[index, 'stored_energy'] = MWh
                df.loc[index, 'revenue_negprices'] = MWh * row.price_negative * -1

for index, row in df.iterrows():        
    energy_storage(df, index, row, MWh=70)  

sns.set(rc={"figure.figsize":(13, 7)})

df['total_storage_revenue'] = df.revenue_negprices + df.revenue_regular
fig, axes = plt.subplots(3, 1, figsize=(13,15), gridspec_kw={'height_ratios': [1, 1,1]})

axes[0] = df.groupby(['year', 'quarter'])[['revenue_negprices','revenue_regular']].sum().plot(kind='bar', stacked=True, ax=axes[0])
axes[0].legend(['Selling storage capacity (Revenue by negative prices)','Re-(selling) stored electricty at favourable prices'])
axes[0].set_title('Revenues by Selling storage capacity & Re-(selling) stored electricity at favourable prices')
axes[0].set_ylabel('Revenue in AUD')

df['total_storage_revenue'] = df.revenue_negprices + df.revenue_regular
axes[1] = df.groupby(['year'])[['revenue_negprices','revenue_regular']].sum().plot(kind='bar', stacked=True, ax=axes[1])
axes[1].legend(['Selling storage capacity (Revenue by negative prices)','Re-(selling) stored electricty at favourable prices'])
axes[1].set_title('Revenues by Selling storage capacity & Re-(selling) stored electricity at favourable prices')
#ax = plt.gca() # get axis
axes[1].set_xticks([0,1,2,3,4,5])
axes[1].set_xticklabels(['2015', '2016', '2017', '2018', '2019', '2020 YTD'], rotation=45) 
axes[1].set_ylabel('Revenue in AUD')

rev_yearly = df.groupby(['year'])[['revenue_negprices','revenue_regular']].sum()
rev_yearly['total'] = rev_yearly.revenue_negprices + rev_yearly.revenue_regular
axes[2] = sns.boxplot(data=rev_yearly, y='total')
axes[2].set_title('Yearly Total Electricity Storage Revenue distribution 2015-2020 YTD')
axes[2].set_ylabel('Total revenue per year in AUD')
plt.show()

print((f"Yearly 25% Percentile: {rev_yearly['total'].describe()[4]:.0f}"))
print((f"Yearly 75% Percentile: {rev_yearly['total'].describe()[6]:.0f}"))
print((f"Yearly Average: {rev_yearly['total'].describe()[1]:.0f}"))

rev_yearly['total'].describe()

from statsmodels.tsa.seasonal import seasonal_decompose
fig2, axes = plt.subplots(2, 1, figsize=(13, 10), constrained_layout=False)
monthly = df.resample('M').mean()
sub1 = sns.barplot(y=monthly.frac_neg_price, x=monthly.index, palette="Blues_d", ax = axes[0])
ticklabels = [monthly.index[int(tick)].strftime('%Y-%m') for tick in sub1.get_xticks()]
sub1.set_xticklabels(ticklabels,rotation=90)
sub1.set_title('monthly mean')
decomp = seasonal_decompose(monthly.frac_neg_price,period=12)
sub2 = decomp.trend.plot(ax = axes[1])
sub2.set_title('isolated monthly trend component pointing upwards')
plt.suptitle('Monthly average of daily fractions, in which transactions traded at negative prices')
plt.subplots_adjust(hspace = 0.5)
plt.show()