<jupyter_start><jupyter_text>Import packages<jupyter_code>import pickle
import datetime as dt
import time
import pandas as pd
import numpy as np
from scipy import stats
import holidays
import pandas_market_calendars as mcal   
from openpyxl import load_workbook
import matplotlib.pyplot as plt
import itertools
from sklearn import linear_model
from sklearn import metrics
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
import seaborn as sns
from sklearn.preprocessing import OneHotEncoder
import sys<jupyter_output><empty_output><jupyter_text>Set display options for dataframes in the notebook<jupyter_code>pd.set_option('display.max_columns', None)  
pd.set_option('display.expand_frame_repr', False)
pd.set_option('max_colwidth', None)<jupyter_output><empty_output><jupyter_text>Read data from file<jupyter_code>df = pd.read_excel ('raw_data.xlsx', index_col = 'dates')<jupyter_output><empty_output><jupyter_text>Separate x and y variables<jupyter_code>y_df = df[['P_wti_oil']].copy()
y_df.dropna(inplace = True)
y_df_pct_chg = y_df.pct_change()
y_df_pct_chg.dropna(inplace = True)
x_df_raw = df.drop(['P_wti_oil'], 1)<jupyter_output><empty_output><jupyter_text>Remove columns with only 0s or NaNs<jupyter_code>col_names = x_df_raw.columns

for name in col_names:
    
    if any(x_df_raw[name].values>0) == False and any(x_df_raw[name].values<0) == False:
        x_df_raw.drop(name, 1, inplace = True)
        print(['column ' + name + ' was removed as there were only non-zero or Nan values.'])<jupyter_output><empty_output><jupyter_text>Separate daily and weekly variables, put weekly variables in a dictionary containing two dataframes: one for positions reports and one for inventory report<jupyter_code>weekly_data = ['pos', 'inventory']
dict_weekly_vars, dict_weekly_vars_pct_chg, dict_weekly_vars_1diff = {}, {}, {}
weekly_cols = list()

for name in weekly_data:
    temp_col = [col for col in x_df_raw.columns if name in col]
    dict_weekly_vars[name] = x_df_raw[temp_col].copy()
    dict_weekly_vars[name].dropna(inplace = True)
    dict_weekly_vars_pct_chg[name] = dict_weekly_vars[name].pct_change()
    dict_weekly_vars_1diff[name] = dict_weekly_vars[name].diff()
    weekly_cols.append(temp_col)

weekly_cols = list(itertools.chain.from_iterable(weekly_cols))
daily_cols = [col for col in x_df_raw.columns if col not in weekly_cols]
x_df_raw_daily = x_df_raw[daily_cols]<jupyter_output><empty_output><jupyter_text>Calculate Net position for positions variables (long minus short) as that's what is useful to know<jupyter_code>df_pos_raw = dict_weekly_vars['pos']

df_pos_raw_net, df_pos_pct_chg_net, df_pos_diff_net = pd.DataFrame(), pd.DataFrame(), pd.DataFrame()

pos_cols = df_pos_raw.columns
pos_cols = [x[:-1] for x in pos_cols]
unique_cat = np.unique(pos_cols)
for i in range(0,4):
    df_pos_raw_net[pos_cols[i] + '_net'] = df_pos_raw.iloc[:,i] - df_pos_raw.iloc[:, i + 4]
    df_pos_pct_chg_net[pos_cols[i] + '_net'] = df_pos_raw_net[pos_cols[i] + '_net'].pct_change()
    df_pos_diff_net[pos_cols[i] + '_net'] = df_pos_raw_net[pos_cols[i] + '_net'].diff()
    is_net_long = df_pos_raw_net[pos_cols[i] + '_net'] > 0
    df_pos_raw_net[pos_cols[i] + '_is_long'] = np.full((len(df_pos_raw_net), 1), 0)
    idx = np.where(is_net_long == True)
    idx = idx[0]
    col_nbr =  df_pos_raw_net.columns.get_loc(pos_cols[i] + '_is_long')
    df_pos_raw_net.iloc[idx, col_nbr] = 1
    
# only keep net positions variables (remove short and long)
dict_weekly_vars['pos'] = df_pos_raw_net
dict_weekly_vars_pct_chg['pos'] = df_pos_pct_chg_net
dict_weekly_vars_diff['pos'] = df_pos_diff_net<jupyter_output><empty_output><jupyter_text>Shift x values by one day from y values as only closing data is available (and therefore usable to make prediction). It's priced in after closing time (the day after).<jupyter_code>x_df_daily = x_df_raw_daily.iloc[0:-1].values
new_idx = x_df_raw_daily.index[1:len(x_df_raw_daily)]
x_df_daily = pd.DataFrame(x_df_daily, index=new_idx, columns=list(x_df_raw_daily.columns))<jupyter_output><empty_output><jupyter_text>Replace missing values by previous value for variables (in our dataset missing values are not due to an error but to the fact that the data point is not available   so we use the previous value instead as that is the piece of information available to the market at the time of the missing value)<jupyter_code>x_df_daily = x_df_daily.fillna(method='ffill')<jupyter_output><empty_output><jupyter_text>Remove rows with no y data<jupyter_code>idx = np.isin(x_df_daily.index, y_df.index)
x_df_daily = x_df_daily.iloc[idx,:]<jupyter_output><empty_output><jupyter_text>Transform data into % change and first difference to make it stationary<jupyter_code>x_df_daily_pct_chg = x_df_daily.pct_change()
x_df_daily_pct_chg.drop(x_df_daily_pct_chg.index[0], inplace = True)
x_df_daily_diff = x_df_daily.diff()
x_df_daily_diff.drop(x_df_daily_diff.index[0], inplace = True)<jupyter_output><empty_output><jupyter_text>Change dates for weekly data to the day they become known to market participants + +5 days for DOE + +3 days for CFCD<jupyter_code>dict_weekly_vars_pct_chg['pos'].index = dict_weekly_vars_pct_chg['pos'].index + pd.DateOffset(days=3)
dict_weekly_vars['pos'].index = dict_weekly_vars['pos'].index + pd.DateOffset(days=3)
dict_weekly_vars_pct_chg['inventory'].index = dict_weekly_vars_pct_chg['inventory'].index + pd.DateOffset(days=5) 
dict_weekly_vars['inventory'].index = dict_weekly_vars['inventory'].index + pd.DateOffset(days=3)<jupyter_output><empty_output><jupyter_text>If the usual weekly variables publication day is a bank holiday, set the publication day to the next market open day<jupyter_code>weekly_vars = ['pos', 'inventory']
for each_var in weekly_vars:
    is_close_week_var =  dict_weekly_vars_pct_chg[each_var].index[~np.isin(dict_weekly_vars_pct_chg[each_var].index, y_df.index)]  
    for each_date in is_close_week_var:
        diff = y_df.index - each_date
        min = diff[(diff > pd.to_timedelta(0))].min()
        each_date_loc = dict_weekly_vars_pct_chg[each_var].index.get_loc(each_date)
        dict_weekly_vars_pct_chg[each_var].index.values[each_date_loc] = each_date + min
    dict_weekly_vars[each_var].index  = dict_weekly_vars_pct_chg[each_var].index<jupyter_output><empty_output><jupyter_text>Handling outliers Outliers visualisation<jupyter_code>x_df_daily_pct_chg_float = x_df_daily_pct_chg.select_dtypes('float64')

sns.set(color_codes=True)
names = list(x_df_daily_pct_chg_float.columns)
f, axes = plt.subplots(len(names), figsize=(10,100))
f2, axes2 = plt.subplots(len(names), figsize=(10,100))
i = 0

for name in names:
    axes[i].scatter(x=x_df_daily_pct_chg_float[name].index, y=x_df_daily_pct_chg_float[name].values)
    axes[i].label_outer()
    axes[i].set_xlabel('years')
    axes[i].title.set_text(name)
    
    sns.boxplot(x=x_df_daily_pct_chg_float[name], ax=axes2[i])
    axes2[i].label_outer()
    axes2[i].set_xlabel('values')
    axes2[i].title.set_text(name)        
    i = i + 1
    
f.tight_layout(pad = 5.0)
f2.tight_layout(pad = 5.0)<jupyter_output><empty_output><jupyter_text>Check for outliers programmaticaly - return dataframe containing outlier analysis<jupyter_code>dict_outliers = {}
# set a threshold for the number of standard deviations above which the value is considered as an outlier
THRESHOLD = 6

for y in df.columns:
    # using 365 days rolling average to calculate z_score as data isn't stationary
    roll_mean = df[y].rolling(365, center=True, min_periods=1).mean()
    roll_std = df[y].rolling(365, center=True, min_periods=1).std()
    dict_outliers[y] = df[[y]].copy()
    dict_outliers[y]['roll_mean'] = roll_mean.values
    dict_outliers[y]['roll_std'] = roll_std.values
    dict_outliers[y]['z_scores'] = abs(dict_outliers[y][y].values - dict_outliers[y]['roll_mean'])/dict_outliers[y]['roll_std']
    dict_outliers[y]['is_outlier'] = False
    
    idx2 = np.where(dict_outliers[y]['z_scores'].values > THRESHOLD)
    idx2 = idx2[0]
    dict_outliers[y].loc[dict_outliers[y].index[idx2], 'is_outlier'] = True<jupyter_output><empty_output><jupyter_text>Remove rows containing outliers<jupyter_code>all_keys = list(dict_outliers.keys())
rows_to_remove = dict_outliers[all_keys[0]]['is_outlier'].index[dict_outliers[all_keys[0]]['is_outlier'] == True]

for key in all_keys[1:]:
    rows_to_remove = rows_to_remove.append(dict_outliers[key]['is_outlier'].index[dict_outliers[key]['is_outlier'] == True])

rows_to_remove.drop_duplicates(keep='first')

# remove outliers in x and y vars in pct as well as raw x and y vars (original unit)

x_df_daily_pct_chg.drop(x_df_daily_pct_chg[x_df_daily_pct_chg.index.isin(rows_to_remove)].index, inplace = True)
x_df_daily.drop(x_df_daily_pct_chg[x_df_daily_pct_chg.index.isin(rows_to_remove)].index, inplace = True)
y_df_pct_chg.drop(y_df_pct_chg[~y_df_pct_chg.index.isin(x_df_daily_pct_chg.index)].index, inplace = True)
y_df.drop(x_df_daily_pct_chg[x_df_daily_pct_chg.index.isin(rows_to_remove)].index, inplace = True)<jupyter_output><empty_output><jupyter_text>Check for stationarity - unit root test<jupyter_code>x_df_daily_pct_chg_adf = x_df_daily_pct_chg.copy()
x_df_daily_pct_chg_adf = x_df_daily_pct_chg_adf.select_dtypes('float64')   
for name in x_df_daily_pct_chg.columns:
    x_serie_daily_pct_chg = x_df_daily_pct_chg_adf[name].dropna()
    try:
        df_test = adfuller(x_serie_daily_pct_chg, autolag='AIC')
        df_result = pd.Series(df_test[0:4], index=['t-stat','p-value','nbr-lags','nbr-obs'])
        for key,value in df_test[4].items():
            df_result['Critical Value (%s)'%key] = value
            if df_result['t-stat'] > df_test[4]['1%']:
                print(['Variable ' + name + ' data is not stationary: Test Statistic = ' + str(dftest[0]) + '; Critical Value 1% = ' + str(dftest[4]['1%'])])
    except:
        the_type, the_value, the_traceback = sys.exc_info()
        print(['Error trying DFT with ' + name + ' ' + str(the_type) + '. ' + str(the_value)])
        continue<jupyter_output><empty_output><jupyter_text>Create lag variables and moving average variables Define function to create lag and mvg average variables:<jupyter_code>def create_lag_mov_avg_var(x_df, get_mvg_avg = False):  
    
    """
    Returns a dataframe containing lagged variables and moving average variables:
    - lagged variables from 1 to 4 periods (days if input dataframe is in days, week if input dataframe is in weeks)
    - moving averages variables over 2 to 10 days (only if get_mvg_avg is set to True)
    """
    
    NUMBER_LAGS = 4
    
    dict_lag_mvavg_vars = {}
    
    for lag in range(1, NUMBER_LAGS + 1):
        temp_val = x_df.iloc[0:-lag].values
        temp_index = x_df.index[lag:len(x_df)]
        dict_lag_mvavg_vars['x_df_daily_lag_' + str(lag)] = pd.DataFrame(temp_val, index=temp_index, columns=list(x_df.columns))
        dict_lag_mvavg_vars['x_df_daily_lag_' + str(lag)] = dict_lag_mvavg_vars['x_df_daily_lag_' + str(lag)].add_suffix('_lag_' + str(lag) + '_p')
        
    # get_mvg_avg is a boolean: true for daily variables for which we want to calculate moving averages 
    # and False for weekly variables for which we don't calculate moving average
    
    if get_mvg_avg:
        # create 2-10 days moving averages
        
        ROLLING_DAYS = 10
        
        for rolling_day in range(2, ROLLING_DAYS + 1):
            dict_lag_mvavg_vars['x_df_rolling_' + str(rolling_day)] = x_df.rolling(rolling_day, min_periods=1).mean()
            dict_lag_mvavg_vars['x_df_rolling_' + str(rolling_day)] = dict_lag_mvavg_vars['x_df_rolling_' + str(rolling_day)].add_suffix('_rolling_avg_' + str(rolling_day) + '_p')
    
    x_df_lag = pd.DataFrame(index = x_df.index)
    
    for key in dict_lag_mvavg_vars.keys():
        x_df_lag = pd.merge(x_df_lag, dict_lag_mvavg_vars[key], how='left', left_index=True, right_index=True)    
        
    return x_df_lag<jupyter_output><empty_output><jupyter_text>Run the function create_lag_mov_avg_var() on daily variables dataframes (variables in pct: x_df_daily_pct_chg and raw values: x_df_daily) and weekly variables (in pct: dict_weekly_vars_pct_chg and raw values: dict_weekly_vars)<jupyter_code>x_df_lag_pct_chg = create_lag_mov_avg_var(x_df_daily_pct_chg, get_mvg_avg = True)
x_df_lag_daily = create_lag_mov_avg_var(x_df_daily, get_mvg_avg = True)

dict_weekly_lag_pct_chg, dict_weekly_lag = {}, {}

for key in dict_weekly_vars_pct_chg.keys():
    # add 1 to all weekly percentage change as when weekly variables are merged with daily variable, NaN on the days weekly variables are missing will be replaced by 0
    # we need to add 1 to weekly percentage changes so inputted values and actual values are easily differentiable
    dict_weekly_vars_pct_chg[key] = dict_weekly_vars_pct_chg[key] + 1
    dict_weekly_lag_pct_chg[key] = create_lag_mov_avg_var(dict_weekly_vars_pct_chg[key], get_mvg_avg = False)
    dict_weekly_lag[key] = create_lag_mov_avg_var(dict_weekly_vars[key], get_mvg_avg = False)<jupyter_output><empty_output><jupyter_text>Create calendar variables<jupyter_code>def create_cal_vars(x_df, y_df, dict_input):
    
    """
    Returns a dataframe containing the following variables, all hot encoded as dummies:
    - day of the month (0-30) 
    - week of the month (0-4) 
    - day of the week (0-6)
    - month (1-12)
    - year (2000-2020)
    - days before or after bank holidays
    The function also merge daily variables with weekly variables.
    """

    temp_df = pd.DataFrame()
    temp_df['day_of_month'] = x_df.index.day
    temp_df['week_of_month'] = (x_df.index.day - 1) // 7 + 1
    temp_df['day_of_week'] = x_df.index.weekday
    temp_df['month'] = x_df.index.month
    temp_df['year'] = x_df.index.year
    temp_df = temp_df.set_index(x_df.index)
    
    # transform all categorical variables as 1/0 dummies
    for col in temp_df.columns:
        temp_df[col] = temp_df[col].astype('category')
        hot_enc = pd.get_dummies(temp_df, prefix= col + '_')
        hot_enc = hot_enc.set_index(x_df.index)#
        x_df_cal = hot_enc.copy()
        #x_df = pd.merge(x_df, hot_enc, how='left', left_index=True, right_index=True)  
        
    #Find which days are closed bank holidays for NYMEX so we can see which days in the y dataset preceeds a bank holiday
    # Get all business dates (no weekend) between the dataset start date and end date
    
    start_date = x_df.index.min().to_pydatetime()
    end_date = x_df.index.max().to_pydatetime()
    all_business_dates = pd.date_range(start_date, end_date, freq='B').to_pydatetime()
    
    # workout which days in the business dates is a bank holiday (NYMEX was closed so no y data on those days)
    NYMEX_open_days = y_df.index.to_pydatetime()
    is_close = ~np.isin(all_business_dates.astype('datetime64[ns]'), NYMEX_open_days.astype('datetime64[ns]'))
    
    # crete a dummy to identify bank holidays (0 if not a bank holiday and 1 if it's a bank holiday)
    temp_df2 = pd.DataFrame(np.full((len(all_business_dates), 4), 0), index=all_business_dates, \
                            columns=['is_trading_days_before_bank_holiday', 'is_days_before_xmas', \
                                     'is_days_after_xmas', 'is_days_after_NY'])
    idx_close = np.where(is_close==True)
    idx_close = idx_close[0]
    
    # identify the 3 days preceeding a bank holiday with a dummy 
    # (0 if not a day preceding a bank holiday - NYMEX close days - 
    # and 1 if the days are within 3 days before a bank holiday)
    
    DAYS_AROUND_BANK_HOL = 3
    
    for x in idx_close:
        temp_df2['is_trading_days_before_bank_holiday'].iloc[x - DAYS_AROUND_BANK_HOL + 1 : x] = 1
        # identify the 3 days before and after xmas and the 3 days after new year
        year = temp_df2.index[x].year
        
        if temp_df2.index[x].to_pydatetime == dt.datetime(year, 12, 25):
            temp_df2['is_days_before_xmas'].iloc[idx - DAYS_AROUND_BANK_HOL : idx + 1] = 1
            temp_df2['is_days_after_xmas'].iloc[idx : idx + DAYS_AROUND_BANK_HOL + 1] = 1
            
        elif temp_df2.index[x].to_pydatetime == dt.datetime(year, 12, 25):
             temp_df2['is_days_after_NY'].iloc[idx : idx + DAYS_AROUND_BANK_HOL + 1] = 1
                
    x_df_cal = pd.merge(x_df_cal, temp_df2, how='left', left_index = True, right_index = True)
    
    # add x weekly variables to x daily variables to put back all x variables together
    weekly_var = ['pos', 'inventory']
    
    for each_var in weekly_var:  
        x_df = pd.merge(x_df, dict_input[each_var], how='left', left_index=True, right_index=True)
        
    # add a dummy to indicate whether the weekly data is available on a day or not
    x_df['is_pos_data'] = np.full((len(x_df), 1), 0)
    x_df['is_pos_data'][~np.isnan(x_df['pos_oil_prod_lon_net'])] = 1
    x_df['is_inventory_data'] = np.full((len(x_df), 1), 0)
    x_df['is_inventory_data'][~np.isnan(x_df['crude_inventory'])] = 1
    
    # Replace NaNs by 0 in missing weekly daily variables:
    for key in dict_input.keys():
        
        for col in dict_input[key].columns:
            x_df[col] = x_df[col].fillna(0)
            
    return x_df_cal, x_df<jupyter_output><empty_output><jupyter_text>Create calendar variables for x vars dataframe in pct and x vars dataframe with raw values<jupyter_code>x_df_cal_pct_chg, x_df_pct_chg = create_cal_vars(x_df_daily_pct_chg, y_df_pct_chg, dict_weekly_vars_pct_chg)
x_df_cal, x_df_daily = create_cal_vars(x_df_daily, y_df, dict_weekly_vars)<jupyter_output><empty_output><jupyter_text>Create products variables Create a function that outputs product variables<jupyter_code>def create_prod_vars(x_df):
    
    """
    Returns a dataframe containing product variables as follow for price variables:
    - volume & price
    - open interest & price
    - volume, open interest & price
    """

    new_cols = [col.replace('P_', '') for col in x_df.columns if 'P_' in col]

    product_vars = pd.DataFrame(index = x_df.index)

    for new_col in new_cols:
        product_vars[new_col + '_V_P'] = x_df['V_' + new_col] * x_df['P_' + new_col]
        product_vars[new_col + '_OI_P'] = x_df['OI_' + new_col] * x_df['P_' + new_col]
        product_vars[new_col + '_V_OI_P'] = x_df['V_' + new_col] * x_df['OI_' + new_col] * x_df['P_' + new_col]

    return product_vars<jupyter_output><empty_output><jupyter_text>Run the create_prod_vars on x dataframes (pct data and raw data)<jupyter_code>product_vars_pct_chg = create_prod_vars(x_df_daily_pct_chg)
product_vars = create_prod_vars(x_df_daily)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10788314896359145, 'fraction_numerical': 0.00950488669315863, 'mean_word_length': 4.746066029003394, 'pattern_counts': {'":': 0, '<': 94, '<?xml version=': 0, '>': 98, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '266717', 'license': 'MIT', 'n_tokens_mistral': 7160, 'n_tokens_neox': 6942, 'n_words': 1423, 'path': 'main_crude_oil_fcst.ipynb', 'repo_name': 'manishanker/time-series-oil-price-prediction', 'size': '18622'}