Time series analysis on EirGrid intra-day wind-power generation

import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns import statsmodels.api as sm from scipy.stats import ttest_ind # from scipy.stats import f_oneway # from arch.unitroot import VarianceRatio # from sklearn.model_selection import train_test_split sns.set_theme(style="darkgrid")

# loading the data grid_data = pd.read_csv('dataset/EirGridSystemDemand2014.csv')

# shape grid_data.shape

# first 5 rows grid_data.head()

Dateobject

Timeobject

0

1/1/2014

0:00

1

1/1/2014

0:15

2

1/1/2014

0:30

3

1/1/2014

0:45

4

1/1/2014

1:00

# last 5 rows grid_data.tail()

Dateobject

Timeobject

35035

12/31/2014

22:45

35036

12/31/2014

23:00

35037

12/31/2014

23:15

35038

12/31/2014

23:30

35039

12/31/2014

23:45

# structure information of the data grid_data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 35040 entries, 0 to 35039
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   Date    35040 non-null  object 
 1   Time    35040 non-null  object 
 2   Demand  35035 non-null  float64
dtypes: float64(1), object(2)
memory usage: 821.4+ KB

# checking what's wrong with 'Demand' grid_data['Demand'].isnull().sum()

# checking the type of object types grid_data.select_dtypes(include=[np.datetime64])

0

1

2

3

4

5

6

7

8

9

grid_data.select_dtypes(include=[np.number])

Demandfloat64

1665.0 - 4614.0

0

2899

1

2869

2

2827

3

2787

4

2724

5

2687

6

2596

7

2543

8

2483

9

2438

grid_data.describe()

Demandfloat64

count

35035

mean

2943.380534

std

599.651709

min

1665

25%

2426

50%

3014

75%

3378

max

4614

# turning 'Date' column to datetime64 grid_data['Date'] = pd.to_datetime(grid_data['Date'])

grid_data['Date'].dtype

# turning 'Time' column to timedelta64 grid_data['Time'] = pd.to_timedelta(grid_data['Time'] + ':00')

# combining date and time with pop() grid_data['Datetime'] = grid_data.pop('Date') + grid_data.pop('Time') # grid_data['Datetime'] = grid_data['Date'] + grid_data['Time'] grid_data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 35040 entries, 0 to 35039
Data columns (total 2 columns):
 #   Column    Non-Null Count  Dtype         
---  ------    --------------  -----         
 0   Demand    35035 non-null  float64       
 1   Datetime  35040 non-null  datetime64[ns]
dtypes: datetime64[ns](1), float64(1)
memory usage: 547.6 KB

# looking for outliers plt.figure(figsize=(12, 6)) sns.boxplot(grid_data['Demand'])

/shared-libs/python3.7/py/lib/python3.7/site-packages/seaborn/_decorators.py:43: FutureWarning: Pass the following variable as a keyword arg: x. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  FutureWarning

# relationship between 'Demand' and 'Datetime' plt.figure(figsize=(12, 6)) sns.scatterplot(grid_data['Datetime'], grid_data['Demand'], label='Demand vs Datetime').set_title('Relationship between Datetime and Demand')

/shared-libs/python3.7/py/lib/python3.7/site-packages/seaborn/_decorators.py:43: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  FutureWarning

# grid_data_og = pd.read_excel('dataset/System-Data-Qtr-Hourly-2014-2015.xlsx') # to csv # grid_data_og.to_csv('dataset/System-Data-Qtr-Hourly-2014-2015.csv', index=False)

# read the original data grid_data_og = pd.read_csv('dataset/System-Data-Qtr-Hourly-2014-2015.csv') grid_data_og.columns

/shared-libs/python3.7/py-core/lib/python3.7/site-packages/IPython/core/interactiveshell.py:3524: DtypeWarning: Columns (12) have mixed types.Specify dtype option on import or set low_memory=False.
  exec(code_obj, self.user_global_ns, self.user_ns)

# drop everything except 'Datetime' and 'Demand' grid_data_og = grid_data_og.drop(grid_data_og.columns.difference(['DateTime', 'IE Demand']), 1) # grid_data_og.drop(grid_data_og.columns.difference(['DateTime', 'IE Demand']), 1, inplace=True) grid_data_og.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 70080 entries, 0 to 70079
Data columns (total 2 columns):
 #   Column     Non-Null Count  Dtype  
---  ------     --------------  -----  
 0   DateTime   70080 non-null  object 
 1   IE Demand  70080 non-null  float64
dtypes: float64(1), object(1)
memory usage: 1.1+ MB
/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:2: FutureWarning: In a future version of pandas all arguments of DataFrame.drop except for the argument 'labels' will be keyword-only

# setting date as datetime64 grid_data_og['DateTime'] = pd.to_datetime(grid_data_og['DateTime']) # select only the data from 2014 grid_data_og = grid_data_og[grid_data_og['DateTime'].dt.year == 2014] grid_data_og.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 35040 entries, 0 to 35039
Data columns (total 2 columns):
 #   Column     Non-Null Count  Dtype         
---  ------     --------------  -----         
 0   DateTime   35040 non-null  datetime64[ns]
 1   IE Demand  35040 non-null  float64       
dtypes: datetime64[ns](1), float64(1)
memory usage: 821.2 KB

# range of dates min_date = grid_data_og['DateTime'].min() max_date = grid_data_og['DateTime'].max() # compute the number of days between the first and last date with a frequency of 15 min range_dates = pd.date_range(start=min_date, end=max_date, freq='15min')

# difference -> look for missing dates in the original data # probably because of the daylight saving time print(range_dates.difference(grid_data_og['DateTime']))

DatetimeIndex(['2014-03-30 01:00:00', '2014-03-30 01:15:00',
               '2014-03-30 01:30:00', '2014-03-30 01:45:00'],
              dtype='datetime64[ns]', freq=None)

# displaying the missing values in the new dataset grid_data[grid_data['Demand'].isna() == True]

Demandfloat64

Datetimedatetime64[ns]

8452

nan

2014-03-30T01:00:00.000000

8453

nan

2014-03-30T01:15:00.000000

8454

nan

2014-03-30T01:30:00.000000

8455

nan

2014-03-30T01:45:00.000000

21657

nan

2014-08-14T14:15:00.000000

# 2014-08-14 14:15:00 of the original data grid_data_og[grid_data_og['DateTime'] == '2014-08-14 14:15:00']

DateTimedatetime64[ns]

IE Demandfloat64

21653

2014-08-14T14:15:00.000000

3270.32

# interpolating missing values grid_data['Demand'] = grid_data['Demand'].interpolate(method='linear')

grid_data.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 35040 entries, 0 to 35039
Data columns (total 2 columns):
 #   Column    Non-Null Count  Dtype         
---  ------    --------------  -----         
 0   Demand    35040 non-null  float64       
 1   Datetime  35040 non-null  datetime64[ns]
dtypes: datetime64[ns](1), float64(1)
memory usage: 547.6 KB

# plot the time series plt.figure(figsize=(15,10)) sns.lineplot(x='Datetime', y='Demand', data=grid_data, label='Demand').set_title('Energy demand overtime') plt.show()

# estimate autocorrelation coefficient for 10 days ten_days_lags_samples = 10 * 24 * 4 fig, ax = plt.subplots(figsize=(15,10)) plt.xlabel('10 days lag') plt.ylabel('Autocorrelation') sm.graphics.tsa.plot_acf(grid_data['Demand'], lags=ten_days_lags_samples, ax=ax, title='Autocorrelation of Energy demand') plt.show()

# 3 # time of year variable that varies from 0 to 1 time_of_year = grid_data['Datetime'].dt.dayofyear / 365 time_of_year

# how demand varies over the year plt.figure(figsize=(15,10)) plt.xlabel("Time of the year(0-1)") plt.ylabel("Demand variation") sns.lineplot(x=time_of_year, y=grid_data['Demand'], data=grid_data, label='Demand variation').set_title('Demand variation overtime') plt.show()

# 4 # average demand per month # avg_demand_per_month = grid_data.resample('M', on='Datetime').mean() avg_demand_per_month = grid_data.groupby(grid_data['Datetime'].dt.month).mean() avg_demand_per_month

Demandfloat64

2681.5524193548385 - 3247.4103422619046

1

3226.265793

2

3247.410342

3

3093.453629

4

2844.273264

5

2767.46203

6

2695.365278

7

2695.177083

8

2681.552419

9

2806.761111

10

2933.15121

# 4 # months = grid_data['Datetime'].dt.month # divide months in 12 categories # months_categories = pd.cut(months, 12, labels=False) twelve_months = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'] # plot average demand per month over months plt.figure(figsize=(15,10)) # bar plot fig = sns.barplot(x=twelve_months, y=avg_demand_per_month['Demand'], data=avg_demand_per_month).set_title('Average Demand per month') plt.xlabel("Months") plt.ylabel("Average Demand") plt.show()

# 5 # average demand each hour of the day avg_demand_per_hour = grid_data.groupby(grid_data['Datetime'].dt.hour).mean() # avg_demand_per_hour = grid_data.resample('H', on='Datetime').mean() # avg_demand_per_hour # # plot average demand per hour over hours plt.figure(figsize=(15,10)) # bar plot sns.barplot(x=avg_demand_per_hour.index, y=avg_demand_per_hour['Demand'], data=avg_demand_per_hour).set_title('Average Demand per hour') plt.xlabel("Hours") plt.ylabel("Average Demand") plt.show()

# 6 # average demand for each day of the week avg_demand_per_day_of_week = grid_data.groupby(grid_data['Datetime'].dt.dayofweek).mean() # avg_demand_per_week = grid_data.resample('D', on='Datetime').mean() avg_demand_per_day_of_week

Demandfloat64

0

2966.307292

1

3052.404647

2

3059.144458

3

3063.700521

4

3031.41847

5

2783.925681

6

2644.316707

# 6 days_of_the_week = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday'] # plot average demand per day of the week plt.figure(figsize=(15,10)) # bar plot sns.barplot(x=days_of_the_week, y=avg_demand_per_day_of_week['Demand'], data=avg_demand_per_day_of_week).set_title('Average Demand per day of the week') plt.xlabel("Days of the week") plt.ylabel("Average demand") plt.show()

# 7 # days_of_the_week_in_year = grid_data['Datetime'].dt.dayofweek # print(days_of_the_week_in_year) # avg_demand_per_hour_of_day = [] # for d in days_of_the_week_in_year: def avg_daily_demand(d): avg_demand_per_hour_of_day = grid_data[grid_data['Datetime'].dt.dayofweek == d].groupby(grid_data['Datetime'].dt.hour).mean() return avg_demand_per_hour_of_day

# print(avg_daily_demand(0)['Demand'].values)

# average daily demand for seven days of the week plt.figure(figsize=(15,10)) for d in range(0, 7): sns.lineplot(x=np.arange(0,24), y=avg_daily_demand(d)['Demand'].values, label=d).set_title('Average Demand per hour in 7 days') plt.xlabel("Hours") plt.ylabel("Average Demand") plt.show()

# 8 # weekdays weekdays = grid_data[grid_data['Datetime'].dt.dayofweek <= 4] weekdays

Demandfloat64

1716.0 - 4614.0

Datetimedatetime64[ns]

2014-01-01 00:00:00 - 2014-12-31 23:45:00

0

2899

2014-01-01T00:00:00.000000

1

2869

2014-01-01T00:15:00.000000

2

2827

2014-01-01T00:30:00.000000

3

2787

2014-01-01T00:45:00.000000

4

2724

2014-01-01T01:00:00.000000

5

2687

2014-01-01T01:15:00.000000

6

2596

2014-01-01T01:30:00.000000

7

2543

2014-01-01T01:45:00.000000

8

2483

2014-01-01T02:00:00.000000

9

2438

2014-01-01T02:15:00.000000

# weekends weekends = grid_data[grid_data['Datetime'].dt.dayofweek > 4] weekends

Demandfloat64

1665.0 - 4152.0

Datetimedatetime64[ns]

2014-01-04 00:00:00 - 2014-12-28 23:45:00

288

3051

2014-01-04T00:00:00.000000

289

2985

2014-01-04T00:15:00.000000

290

2897

2014-01-04T00:30:00.000000

291

2815

2014-01-04T00:45:00.000000

292

2743

2014-01-04T01:00:00.000000

293

2688

2014-01-04T01:15:00.000000

294

2606

2014-01-04T01:30:00.000000

295

2543

2014-01-04T01:45:00.000000

296

2496

2014-01-04T02:00:00.000000

297

2426

2014-01-04T02:15:00.000000

# statistical hypothesis test between weekdays and weekends t_stat, p_value = ttest_ind(weekdays['Demand'], weekends['Demand']) print(t_stat, p_value)

46.54684714308395 0.0

# t_stat, p_value = f_oneway(weekdays['Demand'], weekends['Demand']) # print(t_stat, p_value)

# 9 # split the data into halves train = grid_data[:int(len(grid_data) / 2)] # evaluation set test = grid_data[int(len(grid_data) / 2):]

# mean absolute error def mean_abs_error(actual, predicted): return np.mean(np.abs(actual - predicted))

# persistence benchmark, forecast horizons persistence_model = [] error = [] fh = 24*4 for n in range(0, fh): persistence_model = pd.Series(test['Demand'].shift(n)) persistence_model.dropna(inplace=True) error.append(mean_abs_error(test['Demand'], persistence_model))

# print(error)

# maximum (worst) error print(max(error), "error(energy)","at ", "horizon", error.index(max(error)))

901.8862751831501 error(energy) at  horizon 48

# minimum (best) error print(error[1], "error(energy)","at ", "horizon", error.index(error[1]))

39.70609053028141 error(energy) at  horizon 1

# plot MAE against time leads plt.figure(figsize=(15,10)) plt.xlabel("Forecast horizons") plt.ylabel("Mean absolute error") sns.lineplot(x=range(0, fh), y=error).set_title('Mean absolute error vs forecast horizons') plt.show()

# MAPE def mean_abs_percentage_error(actual, predicted): return np.mean(np.abs((actual - predicted) / actual)) * 100 p_error = [] p_persistence_model = [] for n in range(0, fh): p_persistence_model = pd.Series(test['Demand'].shift(n)) p_persistence_model.dropna(inplace=True) p_error.append(mean_abs_percentage_error(test['Demand'], p_persistence_model))

# print(p_error)

# maximum (worst) error print(max(p_error), "error(energy)","at", "horizon", p_error.index(max(p_error)))

33.31474995425007 error(energy) at horizon 47

# minimum (best) error # print(min(p_error), "error(energy)","at", "horizon", p_error.index(min(p_error))) print(p_error[1], "error(energy)","at ", "horizon", p_error.index(p_error[1]))

1.4014870843965397 error(energy) at  horizon 1

# plot MAPE against time leads plt.figure(figsize=(15,10)) plt.xlabel("Forecast horizons") plt.ylabel("Mean absolute percentage error") sns.lineplot(x=range(0, fh), y=p_error, ).set_title('Mean absolute percentage error vs forecast horizons')