SARIMA Model for Bias Prediction

Here we study and apply ARIMA and SARIMA models, for external model bias prediction.
During the analysis bias can appear as vies, portuguese word for bias.

We present here the equation for SARIMA models:

$y_t = \color{orange}{y_t - d \cdot y_{t-1}} = c + \color{blue}{\sum_{n=1}^{p} \alpha_n y_{t-n}} + \color{red}{\sum_{n=1}^{q} \theta_n \epsilon_{t-n}} + \color{green}{\sum_{n=1}^{P} \phi_n y_{t-sn} + \sum_{n=1}^{Q} \eta_n \epsilon_{t-sn}} + \epsilon_t$

• Blue: AR (Auto Regressive) term (p)
• Red: MA (Moving Average) term (q)
• Orange: Derivative term (d)
• Green: Sazonality AR and MA terms (P, D, Q, s)

# Import libraries
import os
import pandas as pd
import numpy as np
from datetime import date, timedelta
import seaborn as sns
import matplotlib.pyplot as plt
import re

# Modeling libraries
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.seasonal import seasonal_decompose

# Read data function
def read_wind_data(dir_name='Datathons_pem_vento', usina=None):
    """
    Function to read all files containing wind data, returns a dataframe containing all data.
    It considers a fixed frequency of 30 min between each row.
    The folder must contain only wind data, in txt format.
    Each file must follow the pattern UX_Ven_Prev.txt and UX_Ven_Verif.txt, where X is the power plant number.

    Args:
        dir_name (str, optional): string of the folder name, which contains all data
        usina (str, or None): string in UX format, where X is the power plant number.
            It can be None, to read all data inside the folder.

    Returns:
        df: pandas dataframe containing the information of the read data.
            Each column is the name of a file and the row index are ordered timestamps.
    """
    path = os.path.join(os.getcwd(), dir_name)
    df = pd.DataFrame()
    for root, dirs, files in os.walk(path):
        for file in files:
            if usina:
                if re.search(usina + '_Ven_Prev.txt', file) or re.search(usina + '_Ven_Verif.txt', file):
                    file_path = os.path.join(root, file)
                    df_temp = pd.read_csv(file_path, header=None, sep=';')
                    df_temp[0] = pd.to_datetime(df_temp[0], format='%Y%m%d')
                    df_temp.set_index(0, drop=True, inplace=True)
                    init_date = df_temp.index[0]
                    last_date = df_temp.index[-1] + timedelta(days=1) - timedelta(minutes=30)
                    df_temp = df_temp.stack(dropna=False)
                    df_temp.index = pd.date_range(init_date, last_date, freq='30T')
                    df = pd.concat([df, pd.DataFrame(df_temp, columns=[file])], axis=1)
            else:
                file_path = os.path.join(root, file)
                df_temp = pd.read_csv(file_path, header=None, sep=';')
                df_temp[0] = pd.to_datetime(df_temp[0], format='%Y%m%d')
                df_temp.set_index(0, drop=True, inplace=True)
                init_date = df_temp.index[0]
                last_date = df_temp.index[-1] + timedelta(days=1) - timedelta(minutes=30)
                df_temp = df_temp.stack(dropna=False)
                df_temp.index = pd.date_range(init_date, last_date, freq='30T')
                df = pd.concat([df, pd.DataFrame(df_temp, columns=[file])], axis=1)
    df = df.reindex(sorted(df.columns), axis=1)
    return df

data = read_wind_data()
data

	U10_Ven_Prev.txt	U10_Ven_Verif.txt	U11_Ven_Prev.txt	U11_Ven_Verif.txt	U12_Ven_Prev.txt	U12_Ven_Verif.txt	U13_Ven_Prev.txt	U13_Ven_Verif.txt	U14_Ven_Prev.txt	U14_Ven_Verif.txt	...	U5_Ven_Prev.txt	U5_Ven_Verif.txt	U6_Ven_Prev.txt	U6_Ven_Verif.txt	U7_Ven_Prev.txt	U7_Ven_Verif.txt	U8_Ven_Prev.txt	U8_Ven_Verif.txt	U9_Ven_Prev.txt	U9_Ven_Verif.txt
2017-01-01 00:00:00	3.8796	3.2332	NaN	1.2080	NaN	1.1268	NaN	1.7056	NaN	1.9205	...	2.9649	3.0499	2.6004	4.4208	4.5385	3.7831	4.7348	3.2173	4.5700	3.2301
2017-01-01 00:30:00	3.8306	3.0635	NaN	1.0023	NaN	0.6833	NaN	1.8107	NaN	2.2813	...	2.9316	2.7315	2.5076	4.4953	4.5385	3.4434	4.5437	3.2036	4.5069	3.1167
2017-01-01 01:00:00	3.7815	2.7986	NaN	0.8058	NaN	0.7950	NaN	3.8786	NaN	2.7579	...	2.8983	2.4952	2.4147	3.8723	4.5385	3.7608	4.3527	3.3689	4.4439	2.7939
2017-01-01 01:30:00	3.7377	3.0123	NaN	1.0192	NaN	1.1210	NaN	3.7506	NaN	2.7623	...	2.9053	2.8445	2.3656	3.5962	4.4579	3.0795	4.1530	3.7987	4.3843	2.7652
2017-01-01 02:00:00	3.6939	2.9413	NaN	1.0742	NaN	1.0895	NaN	2.6184	NaN	2.3249	...	2.9123	2.9139	2.3166	3.5529	4.3773	2.7538	3.9532	3.9098	4.3247	2.8384
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2020-12-31 21:30:00	3.1840	2.8725	2.0309	1.5289	1.0023	1.8960	1.6787	1.2904	1.7821	1.4464	...	2.4182	3.1559	3.1577	4.2607	2.9018	2.7038	3.1121	3.1995	4.2213	2.8216
2020-12-31 22:00:00	3.0701	2.3312	1.9381	1.1669	0.9252	1.6594	1.5701	0.8090	1.6682	1.2822	...	2.2570	3.0932	2.9965	3.9901	3.2137	2.7238	3.4275	3.2473	4.1320	2.6242
2020-12-31 22:30:00	3.0297	2.5118	1.6542	0.7361	0.9042	1.6658	1.4124	0.5454	1.3703	1.2373	...	2.0975	2.7557	2.6828	3.9751	3.3171	3.0499	3.6501	3.6175	3.8726	2.8661
2020-12-31 23:00:00	2.9894	2.9362	1.3703	0.9284	0.8832	1.5034	1.2547	0.9915	1.0724	1.2548	...	1.9381	2.9703	2.3691	3.9582	3.4205	3.1400	3.8726	3.4096	3.6133	3.3746
2020-12-31 23:30:00	2.9386	2.6560	1.0479	0.9755	0.9603	1.5967	1.1846	1.4557	0.9112	1.3022	...	1.8469	2.6815	2.0485	3.5424	3.2453	2.8359	3.7342	3.2759	3.2856	4.3829

70128 rows × 28 columns

ARIMA/SARIMA Model

We will apply both ARIMA and SARIMA models, to predict the bias. Therefore, our time-series will be other prediction models bias. Our study will be applied to only one power plant (U1) and in the end we will see its outcome on other power plants.

df_u1 = data[['U1_Ven_Verif.txt', 'U1_Ven_Prev.txt']]

fig, ax = plt.subplots(figsize=(12,6))
df_u1.plot(ax=ax)
fig.suptitle('Full time series of Power Plant U1')
plt.show()

df_u1.loc[:, 'vies'] = df_u1['U1_Ven_Verif.txt'] - df_u1['U1_Ven_Prev.txt']

/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1597: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.obj[key] = value
/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1676: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self._setitem_single_column(ilocs[0], value, pi)

df_u1

	U1_Ven_Verif.txt	U1_Ven_Prev.txt	vies
2017-01-01 00:00:00	3.2928	2.9088	0.3840
2017-01-01 00:30:00	4.6450	2.7862	1.8588
2017-01-01 01:00:00	4.4030	2.6635	1.7395
2017-01-01 01:30:00	4.4062	2.5356	1.8706
2017-01-01 02:00:00	4.9379	2.4077	2.5302
...	...	...	...
2020-12-31 21:30:00	2.8684	2.5969	0.2715
2020-12-31 22:00:00	2.9722	2.9404	0.0318
2020-12-31 22:30:00	2.8980	2.9930	-0.0950
2020-12-31 23:00:00	2.6601	3.0455	-0.3854
2020-12-31 23:30:00	2.7385	2.8177	-0.0792

70128 rows × 3 columns

# Replacing NAN bias by median
df_u1.loc[:, 'vies'].fillna(df_u1['vies'].median(), inplace=True)

/home/silas/.local/lib/python3.8/site-packages/pandas/core/series.py:4463: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  return super().fillna(

fig, ax = plt.subplots(figsize=(12,6))
df_u1['vies'].plot(ax=ax)
fig.suptitle('Full bias of Power Plant U1')
plt.show()

First Analysis

We start by analyzing data from 1 month, used as training set and prediction of 2 days. We apply window of 1 month for training set, accordingly to the document "Modelo_de_correco_de_vies_da_previso_de_vento_do_ONS", where it is said that, for ONS prediction models, 1 month window showed the best performance. A future work may analyze the impact of the training window over these SARIMA models.

In the first analysis we will apply November/2020 as training set.

df = df_u1[df_u1.index >= '2020-11-01'][['vies']]

fig, ax = plt.subplots(figsize=(12,6))
df.plot(ax=ax)
fig.suptitle('Nov/2020 bias of Power Plant U1')
plt.show()

train = df[df.index < '2020-12-01']
test = df[df.index >= '2020-12-01']

Signal Stationarity Analysis

To guarantee that SARIMA models will work, we need to guarantee the stationarity of the data. This is obtained here, through the mean and standard deviation visualization and the statistical adfuller test.

def adfuller_test(data):
    result = adfuller(data)

    print('Results of Dickey-Fuller Test:')
    print('ADF Statistic: {}'.format(result[0]))
    print('p-value: {}'.format(result[1]))
    print('#Lags used: {}'.format(result[2]))
    print('#Observations used: {}'.format(result[3]))
    print('Critical Values:')
    for key, value in result[4].items():
        print('\t{}: {}'.format(key, value))

# 48 points correspond to 24h, or one day
rolling_mean = train.rolling(window = 48).mean()
rolling_std = train.rolling(window = 48).std()

fig, ax = plt.subplots(figsize=(12,6))
plt.plot(train, color = 'blue', label = 'Original', alpha=0.5)
plt.plot(rolling_mean, color = 'red', label = 'Rolling Mean')
plt.plot(rolling_std, color = 'black', label = 'Rolling Std')
plt.legend(loc = 'best')
plt.title('Rolling Mean & Rolling Standard Deviation')
plt.show()

adfuller_test(train)

Results of Dickey-Fuller Test:
ADF Statistic: -9.50335301139156
p-value: 3.3983433554294535e-16
#Lags used: 16
#Observations used: 1423
Critical Values:
	1%: -3.434953749425684
	5%: -2.8635732362456845
	10%: -2.567852483484606

The visualization shows that the daily mean fluctuates widely, but it does not appear to show a trend. This is confirmed by the adfuller test, which has a p-value of the order of 10e-16.

The signal, however, appears to present a daily pattern, which we will take in account later.

Correlogram and partial correlogram

Our next step is to analyze the ACF and PACF graphics, to identify the ARIMA and SARIMA model hiperparameters.

plot_acf(train, lags=50)
plt.show()

plot_pacf(train, method='ywm', lags=50)
plt.show()

The data autocorrelation clearly shows the presence of sazonality. Before proceeding to this, we will analyze the signals first derivative.

1st derivative

train_diff = train.diff()
train_diff = train_diff.dropna()

# Plot
fig, axes = plt.subplots(2, 1, figsize=(10,5), dpi=100, sharex=True)

# Original
axes[0].plot(train)
axes[0].set_title('Original Series')

# Seasonal
axes[1].plot(train_diff)
axes[1].set_title('1st Derivative')
plt.suptitle('Bias- Nov/2020', fontsize=16)
plt.show()

adfuller_test(train_diff)

Results of Dickey-Fuller Test:
ADF Statistic: -11.32560419250231
p-value: 1.1455618208211097e-20
#Lags used: 24
#Observations used: 1414
Critical Values:
	1%: -3.4349831053858746
	5%: -2.863586191208161
	10%: -2.567859382293453

plot_acf(train_diff, lags=50)
plt.show()

plot_pacf(train_diff, method='ywm', lags=50)
plt.show()

ARIMA (p,d,q) = (0,1,1)

With these results, we apply our first, simpler, model, an ARIMA with derivative.
Below are presented the model parameters as well as some statistics analysis, used to investigate the model quality and reliability.

model = SARIMAX(train, order=(0,1,1))
model = model.fit()

model.summary()

SARIMAX Results

Dep. Variable:	vies	No. Observations:	1440
Model:	SARIMAX(0, 1, 1)	Log Likelihood	-660.664
Date:	Tue, 12 Jul 2022	AIC	1325.329
Time:	16:28:40	BIC	1335.872
Sample:	11-01-2020	HQIC	1329.265
	- 11-30-2020
Covariance Type:	opg

	coef	std err	z	P>|z|	[0.025	0.975]
ma.L1	0.1643	0.020	8.229	0.000	0.125	0.203
sigma2	0.1467	0.003	42.683	0.000	0.140	0.153

Ljung-Box (L1) (Q):	0.02	Jarque-Bera (JB):	612.42
Prob(Q):	0.89	Prob(JB):	0.00
Heteroskedasticity (H):	0.73	Skew:	0.41
Prob(H) (two-sided):	0.00	Kurtosis:	6.09

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

model.plot_diagnostics(figsize=(12,8))
plt.show()

df['forecast'] = model.predict(start='2020-12-01', end='2020-12-10', dynamic=True)
df[['vies','forecast']].plot(figsize=(12,8))
plt.xlim(['2020-11-30', '2020-12-03'])
plt.show()

df

	vies	forecast
2020-11-01 00:00:00	-0.4182	NaN
2020-11-01 00:30:00	0.8832	NaN
2020-11-01 01:00:00	1.4046	NaN
2020-11-01 01:30:00	1.0087	NaN
2020-11-01 02:00:00	0.6152	NaN
...	...	...
2020-12-31 21:30:00	0.2715	NaN
2020-12-31 22:00:00	0.0318	NaN
2020-12-31 22:30:00	-0.0950	NaN
2020-12-31 23:00:00	-0.3854	NaN
2020-12-31 23:30:00	-0.0792	NaN

2928 rows × 2 columns

# Prediction in-sample 1 step ahead
df['forecast'] = model.predict(start='2020-11-01', end='2020-12-01', dynamic=False)
df[['vies','forecast']].plot(figsize=(12,8))
plt.xlim(['2020-11-01', '2020-12-01'])
plt.show()

# Prediction in-sample 1 step ahead
df['forecast'] = model.predict(start='2020-11-01', end='2020-12-01', dynamic=False)
df[['vies','forecast']].plot(figsize=(12,8))
plt.xlim(['2020-11-01', '2020-11-02'])
plt.show()

The in-sample 1 step ahead prediction is used to show that the model is capable of following the data. However, since our data contains sazonality, this simpler model is capable of only predicting the average daily value. We proceed to SARIMA model.

SARIMA Model

For our SARIMA model, we first investigate its data decomposition, in trend, sazonality and residue. The sazonality compromises of 48 points, corresponding to a daily sazonality.

seasonal = seasonal_decompose(df['vies'], model='additive', period=48)
fig = seasonal.plot()
fig.set_size_inches((12, 8))
plt.show()

df = df_u1[df_u1.index >= '2020-11-01'][['vies']]
df.loc[:,'trend'] = df['vies'].diff(48)
df = df.dropna()

train = df[df.index < '2020-12-01']
test = df[df.index >= '2020-12-01']

# Plot
fig, axes = plt.subplots(2, 1, figsize=(10,5), dpi=100, sharex=True)

# Original
axes[0].plot(train['vies'])
axes[0].set_title('Original Series')

# Seasonal
axes[1].plot(train['trend'])
axes[1].set_title('Sazonal Differentiation')
plt.suptitle('Bias Nov/2020', fontsize=16)
plt.show()

Stationarity test

# 48 points correspond to 24h, or one day
rolling_mean = train['trend'].rolling(window = 48).mean()
rolling_std = train['trend'].rolling(window = 48).std()

fig, ax = plt.subplots(figsize=(12,6))
plt.plot(train['trend'], color = 'blue', label = 'Original', alpha=0.5)
plt.plot(rolling_mean, color = 'red', label = 'Rolling Mean')
plt.plot(rolling_std, color = 'black', label = 'Rolling Std')
plt.legend(loc = 'best')
plt.title('Rolling Mean & Rolling Standard Deviation')
plt.show()

adfuller_test(train['trend'])

Results of Dickey-Fuller Test:
ADF Statistic: -7.267452010898051
p-value: 1.6209901372686905e-10
#Lags used: 10
#Observations used: 1381
Critical Values:
	1%: -3.435094023613352
	5%: -2.863635138949194
	10%: -2.567885448250043

Correlogram and partial correlogram

plot_acf(train['trend'], lags=50)
plt.show()

plot_pacf(train['trend'], method='ywm', lags=50)
plt.show()

SARIMA: (p,d,q)(P,D,Q,s) = (1,1,1)(1,1,1,48)

We performed a grid search over some hiperparameters values and here we present the best one achieved.

model_u1 = SARIMAX(train['vies'], order=(1,1,1), seasonal_order=(1,1,1,48))
model_u1 = model_u1.fit()

model_u1.summary()

SARIMAX Results

Dep. Variable:	vies	No. Observations:	1392
Model:	SARIMAX(1, 1, 1)x(1, 1, 1, 48)	Log Likelihood	-661.672
Date:	Tue, 12 Jul 2022	AIC	1333.345
Time:	16:33:56	BIC	1359.358
Sample:	11-02-2020	HQIC	1343.089
	- 11-30-2020
Covariance Type:	opg

	coef	std err	z	P>|z|	[0.025	0.975]
ar.L1	-0.2436	0.177	-1.377	0.169	-0.590	0.103
ma.L1	0.3574	0.170	2.100	0.036	0.024	0.691
ar.S.L48	0.0321	0.029	1.100	0.271	-0.025	0.089
ma.S.L48	-0.9930	0.172	-5.762	0.000	-1.331	-0.655
sigma2	0.1404	0.023	6.201	0.000	0.096	0.185

Ljung-Box (L1) (Q):	0.13	Jarque-Bera (JB):	450.41
Prob(Q):	0.72	Prob(JB):	0.00
Heteroskedasticity (H):	0.81	Skew:	0.38
Prob(H) (two-sided):	0.03	Kurtosis:	5.73

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

model_u1.plot_diagnostics(figsize=(12,8))
plt.show()

df['forecast'] = model_u1.predict(start='2020-12-01', end='2020-12-03', dynamic=True)
df[['vies','forecast']].plot(figsize=(12,8))
plt.xlim(['2020-11-30', '2020-12-03'])
plt.show()

# Previsao in-sample 1 passo a frente
df['forecast'] = model_u1.predict(start='2020-11-10', end='2020-11-12', dynamic=False)
df[['vies','forecast']].plot(figsize=(12,8))
plt.xlim(['2020-11-10', '2020-11-12'])
plt.show()

# NMAPE Calculation:
df.loc[:, 'error'] = np.abs(df['vies'] - df['forecast'])
mean = df[(df.index >= '2020-12-01') & (df.index < '2020-12-03')]['vies'].mean()

nmape = df[(df.index >= '2020-12-01') & (df.index < '2020-12-03')]['error'].sum()
nmape = nmape*100/(96*mean)

print(f'The prediction bias NMAPE error for 2 days window was: {nmape} %')

The prediction bias NMAPE error for 2 days window was: 36.5491606085637 %

Same Model applied to Power Plants 3 and 13

We selected Power Plants 3 and 13 to test this model behavior on other power plants. These two presented worse bias, with higher deviation, among others.

Power Plant 3

df_u3 = data[['U3_Ven_Verif.txt', 'U3_Ven_Prev.txt']]

df_u3.loc[:, 'vies'] = df_u3['U3_Ven_Verif.txt'] - df_u3['U3_Ven_Prev.txt']
# Substituindo vies NAN pela mediana
df_u3.loc[:, 'vies'].fillna(df_u3['vies'].median(), inplace=True)

df = df_u3[df_u3.index >= '2020-11-01'][['vies']]

train = df[df.index < '2020-12-01']
test = df[df.index >= '2020-12-01']

/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1597: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.obj[key] = value
/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1676: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self._setitem_single_column(ilocs[0], value, pi)
/home/silas/.local/lib/python3.8/site-packages/pandas/core/series.py:4463: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  return super().fillna(

model_u3 = SARIMAX(train['vies'], order=(1,1,1), seasonal_order=(1,1,1,48))
model_u3 = model_u3.fit()

df['forecast'] = model_u3.predict(start='2020-12-01', end='2020-12-03', dynamic=True)
df[['vies','forecast']].plot(figsize=(12,8))
plt.xlim(['2020-11-30', '2020-12-03'])
plt.show()

# NMAPE Calculation:
df.loc[:, 'error'] = np.abs(df['vies'] - df['forecast'])
mean = df[(df.index >= '2020-12-01') & (df.index < '2020-12-03')]['vies'].mean()

nmape = df[(df.index >= '2020-12-01') & (df.index < '2020-12-03')]['error'].sum()
nmape = nmape*100/(96*mean)

print(f'The prediction bias NMAPE error for 2 days window was: {nmape} %')

The prediction bias NMAPE error for 2 days window was: 44.63037878080562 %

Power Plant 13

df_u13 = data[['U13_Ven_Verif.txt', 'U13_Ven_Prev.txt']]

df_u13.loc[:, 'vies'] = df_u13['U13_Ven_Verif.txt'] - df_u13['U13_Ven_Prev.txt']
# Substituindo vies NAN pela mediana
df_u13.loc[:, 'vies'].fillna(df_u13['vies'].median(), inplace=True)

df = df_u13[df_u13.index >= '2020-11-01'][['vies']]

train = df[df.index < '2020-12-01']
test = df[df.index >= '2020-12-01']

/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1597: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self.obj[key] = value
/home/silas/.local/lib/python3.8/site-packages/pandas/core/indexing.py:1676: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  self._setitem_single_column(ilocs[0], value, pi)
/home/silas/.local/lib/python3.8/site-packages/pandas/core/series.py:4463: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  return super().fillna(

model_u13 = SARIMAX(train['vies'], order=(1,1,1), seasonal_order=(1,1,1,48))
model_u13 = model_u13.fit()

df['forecast'] = model_u13.predict(start='2020-12-01', end='2020-12-03', dynamic=True)
df[['vies','forecast']].plot(figsize=(12,8))
plt.xlim(['2020-11-30', '2020-12-03'])
plt.show()

# NMAPE Calculation:
df.loc[:, 'error'] = np.abs(df['vies'] - df['forecast'])
mean = df[(df.index >= '2020-12-01') & (df.index < '2020-12-03')]['vies'].mean()

nmape = df[(df.index >= '2020-12-01') & (df.index < '2020-12-03')]['error'].sum()
nmape = nmape*100/(96*mean)

print(f'The prediction bias NMAPE error for 2 days window was: {nmape} %')

The prediction bias NMAPE error for 2 days window was: 101.7233554933708 %

Post work

These results show that a good model for one power plant may not be the best one for others. Each one need to be analyzed individually. Moreover, a full analysis must be conducted, regarding other hyperparameters like window size and different sazonalities.

Our model considered a fixed daily sazonality and addictive SARIMA models. A different approach would be consider a multiplicative SARIMA model. We can also extend to SARIMAX models, using external data like weather prediction as exogenous predictors.