<jupyter_start><jupyter_text>Load carbon intensity data and visualize missing data<jupyter_code># Select the regions we want to examine

ba_list = ['EIA.CISO', 'EIA.ISNE', 'EIA.BPAT','EIA.MISO', 'EIA.NYIS', 'EIA.PJM', 'EIA.SWPP',
'AEC', 'AECI', 'AVA', 'AZPS', 'BANC', 'CHPD', 'CPLE', 'CPLW', 'DOPD', 'DUK', 'EPE', 'ERCO', 'FMPP', 
'FPC', 'FPL', 'GCPD', 'GVL', 'HST','IID', 'IPCO', 'JEA', 'LDWP', 'LGEE', 'NEVP', 'NSB',
'NWMT', 'PACE', 'PACW', 'PGE', 'PNM', 'PSCO', 'PSEI', 'SC', 'SCEG', 'SCL', 'SEC', 
'SOCO', 'SPA', 'SRP', 'TAL', 'TEC', 'TEPC', 'TIDC', 'TPWR', 'TVA', 'WACM', 'WALC', 'WAUW']


ef_year = 2019

# specify the type of emission factor we want to use
ef_type = 'consumption_ef_EGRID_2019'

# choose the units you want to use -either kgCO2/kWh or lbCO2/MWh
units = 'lbCO2/MWh'

# the base units are kg/kWh
unit_conversion = 1
if units == 'lbCO2/MWh':
    unit_conversion = 1 / 0.453592 * 1000

print(f'Year = {ef_year}')
print(f'Number of BAs to load: {len(ba_list)}')

# Load the hourly emission factors
##################################
hourly_ef = load_data.load_hourly_efs(ba_list, ef_year, ef_type)

# rename columns to remove EIA. prefix
hourly_ef.columns = [col.split('.')[-1] for col in hourly_ef.columns]
# update list of regions
ba_list = list(hourly_ef.columns)

# remove negative values
hourly_ef[hourly_ef < 0] = np.NaN

# create a dictionary that matches BA codes to names
ba_name_dict = pd.read_csv('../data/manual/ba_names.csv', index_col='ba_code').to_dict()['ba_name']

print(f'Number of BAs available: {len(ba_list)}')<jupyter_output><empty_output><jupyter_text>Visualize Missing Data<jupyter_code>missingno.matrix(hourly_ef, labels=True)<jupyter_output><empty_output><jupyter_text>Drop BAs with all missing data<jupyter_code>print('BAs with all data missing were dropped:')
print(hourly_ef.isna().sum()[hourly_ef.isna().sum() == 8760])

# drop any columns with all missing values
hourly_ef = hourly_ef.dropna(axis=1, how='all')

#update the BA list
ba_list = list(hourly_ef.columns)

print(f'Number of BAs available: {len(ba_list)}')<jupyter_output><empty_output><jupyter_text>Restore Missing Values from EIA-930Carbonara had filled missing data by filling forward gaps of less than 24 hours. Using the EIA-930 data, we want to identify where the missing source data exists, and use this to restore the missing values in the EF dataset so that we can perform our own interpolation of the missing data.<jupyter_code># get variance statistics for each BA for the Carbonara data(this will be used later)
variance_carbonara = pd.DataFrame()
variance_carbonara['stdev'] = hourly_ef.std(ddof=0)
# load the data from EIA-930 for 2019
columns_to_use = ['Balancing Authority', 'UTC Time at End of Hour',
       'Net Generation (MW) from Coal', 'Net Generation (MW) from Natural Gas',
       'Net Generation (MW) from Nuclear',
       'Net Generation (MW) from All Petroleum Products',
       'Net Generation (MW) from Hydropower and Pumped Storage',
       'Net Generation (MW) from Solar', 'Net Generation (MW) from Wind',
       'Net Generation (MW) from Other Fuel Sources',
       'Net Generation (MW) from Unknown Fuel Sources']

eia_930 = pd.concat([pd.read_csv('../data/downloaded/eia/EIA930_BALANCE_2019_Jan_Jun.csv', usecols=columns_to_use, thousands=',', parse_dates=['UTC Time at End of Hour']),
                     pd.read_csv('../data/downloaded/eia/EIA930_BALANCE_2019_Jul_Dec.csv', usecols=columns_to_use, thousands=',', parse_dates=['UTC Time at End of Hour'])])

#only keep regions that are in the BA list
eia_930 = eia_930[eia_930['Balancing Authority'].isin(ba_list)]

# get a list of unique BAs in this dataset
bas_in_930 = list(eia_930['Balancing Authority'].unique())

# convert from end of hour timestamp to beginning of hour timestamp
eia_930['UTC Time at End of Hour'] = eia_930['UTC Time at End of Hour'] - pd.Timedelta(hours=1)

# localize the timezone as UTC time
eia_930['UTC Time at End of Hour'] = eia_930['UTC Time at End of Hour'].dt.tz_localize('UTC')

# we want to create a datetime column that is in local standard time

# create a column for local time
eia_930['Local Time'] = eia_930['UTC Time at End of Hour']

# convert from UTC time to local time
for ba in bas_in_930:
    # get the time zone
    local_tz = download_data.ba_timezone(ba, 'GMT')

    # populate the local time column, then strip the TZ info from the data
    eia_930.loc[eia_930['Balancing Authority'] == ba, 'Local Time'] = eia_930.loc[eia_930['Balancing Authority'] == ba, 'UTC Time at End of Hour'].dt.tz_convert(local_tz).dt.tz_localize(None)

# drop the UTC column, and set the BA and local time as index
eia_930 = eia_930.drop(columns=['UTC Time at End of Hour'])
eia_930 = eia_930.set_index(['Balancing Authority','Local Time'])

eia_930.head(3)
# find all rows where there is no net generation data
missing_930 = eia_930[eia_930.sum(axis=1) == 0].reset_index()

# get a list of all BAs that have missing net generation data
missing_bas = list(missing_930['Balancing Authority'].unique())

# for each BA, replace the filled values in hourly_ef with missing values where the source data is missing
for ba in missing_bas:
    # get the datetimes that are missing and replace with NaN
    hourly_ef.loc[missing_930[missing_930['Balancing Authority'] == ba]['Local Time'], ba] = np.NaN

# remove negative values
hourly_ef[hourly_ef < 0] = np.NaN

print('BAs with all data missing were dropped:')
print(hourly_ef.isna().sum()[hourly_ef.isna().sum() == 8760])

# drop any columns with all missing values
hourly_ef = hourly_ef.dropna(axis=1, how='all')

#update the BA list
ba_list = list(hourly_ef.columns)

print(f'Number of BAs available: {len(ba_list)}')
missingno.matrix(hourly_ef, labels=True)<jupyter_output><empty_output><jupyter_text>Test the performance of different interpolation methodsWhen filling missing data, we want to be careful to not affect the variance of the data too much, since this could bias our hourly inventory results.We will check several filling methods to see how much this impacts the variance compared to the original data with all the missing values intact.<jupyter_code># calculate variance stats for the restored missing data
variance_missing = pd.DataFrame()
variance_missing['stdev'] = hourly_ef.std(ddof=0)

# calculate the mean absolute percent difference between the raw data and carbonara data
print(f'Mean Absolute percentage difference: {round(abs(((variance_carbonara - variance_missing) / variance_missing)).mean().values[0] *100, 2)}%')
print(f'Mean percentage difference: {round(((variance_carbonara - variance_missing) / variance_missing).mean().values[0] *100, 2)}%')<jupyter_output><empty_output><jupyter_text>Fill-forwardFilling forward all data makes the percent difference worse at 0.34%, which makes sense since we are filling more missing values than Carbonara<jupyter_code>ef_ffill = hourly_ef.copy().fillna(method='ffill')

variance_ffill = pd.DataFrame()
variance_ffill['stdev'] = ef_ffill.std(ddof=0)

print(f'Mean Absolute percentage difference: {round(abs(((variance_ffill - variance_missing) / variance_missing)).mean().values[0] *100, 2)}%')
print(f'Mean percentage difference: {round(((variance_ffill - variance_missing) / variance_missing).mean().values[0] *100, 2)}%')<jupyter_output><empty_output><jupyter_text>Linear InterpolationLinear interpolation is a slight improvement to 0.29%.<jupyter_code>ef_linear = hourly_ef.copy().interpolate(method='linear')

variance_linear = pd.DataFrame()
variance_linear['stdev'] = ef_linear.std(ddof=0)

print(f'Mean Absolute percentage difference: {round(abs(((variance_linear - variance_missing) / variance_missing)).mean().values[0] *100, 2)}%')
print(f'Mean percentage difference: {round(((variance_linear - variance_missing) / variance_missing).mean().values[0] *100, 2)}%')<jupyter_output><empty_output><jupyter_text>Moving AverageUsing a 4-day, centered moving average, which itself is linearly interpolated for larger gaps reduces the percent difference to 0.20%.We tested several windows (24, 36, 48, 72, 96, 120, and 168 hours) and this seemed to perform best<jupyter_code>ef_rolling = hourly_ef.copy().fillna(hourly_ef.rolling(window=96, center=True, min_periods=24, axis=0).mean().interpolate(method='linear'))

variance_rolling = pd.DataFrame()
variance_rolling['stdev'] = ef_rolling.std(ddof=0)

print(f'Mean Absolute percentage difference: {round(abs(((variance_rolling - variance_missing) / variance_missing)).mean().values[0] *100, 2)}%')
print(f'Mean percentage difference: {round(((variance_rolling - variance_missing) / variance_missing).mean().values[0] *100, 2)}%')<jupyter_output><empty_output><jupyter_text>Moving average plus hourly fixed effectCould adding an hourly fixed effect to the moving average improve this further?For each BA that has missing values, calculate the moving average, and also construct an hourly fixed effect model for the whole year. Then add the hourly fixed effect to the moving average.<jupyter_code>ef_hourly_MA = hourly_ef.copy()

# for each BA that has any missing values
for ba in list(ef_hourly_MA.columns):
    if ef_hourly_MA[ba].isna().any():
        # create a df with only that BA
        hourly_model = ef_hourly_MA.copy()[[ba]]
        # ad a column for the hour
        hourly_model['hour'] = hourly_model.index.hour
        # calculate a model
        model = ols(f'{ba} ~ C(hour)', data=hourly_model).fit()
        # get a list of the coefficients
        params = model.params.values
        # replace the first value (the intercept) with 0
        params[0] = 0.0
        # add these values to the df
        hourly_model['fixed_effect'] = np.tile(params, 365)
        # calculate the moving average, and linearly interpolate any gaps in this data
        hourly_model['MA'] = ef_hourly_MA[ba].rolling(window=96, center=True, min_periods=24, axis=0).mean().interpolate(method='linear')
        hourly_model['filled'] = hourly_model[ba].fillna(hourly_model['MA'] + hourly_model['fixed_effect'])

        # fill the data with the missing values
        ef_hourly_MA[ba] = ef_hourly_MA[ba].fillna(hourly_model['filled'])

variance_hourly_MA = pd.DataFrame()
variance_hourly_MA['stdev'] = ef_hourly_MA.std(ddof=0)

print(f'Mean Absolute percentage difference: {round(abs(((variance_hourly_MA - variance_missing) / variance_missing)).mean().values[0] *100, 2)}%')
print(f'Mean percentage difference: {round(((variance_hourly_MA - variance_missing) / variance_missing).mean().values[0] *100, 2)}%')<jupyter_output><empty_output><jupyter_text>monthly TOD averageFill with the monthly time of day average<jupyter_code>ef_mh = hourly_ef.copy()

# calculate the month-hour average
mh_average = ef_mh.copy()
mh_average['month'] = mh_average.index.month
mh_average['hour'] = mh_average.index.hour
mh_average = mh_average.groupby(['month','hour']).mean().reset_index()

mh_hold = pd.DataFrame(index=ef_mh.index)
mh_hold['month'] = mh_hold.index.month
mh_hold['hour'] = mh_hold.index.hour

# merge month-hourly
mh_hold = mh_hold.merge(mh_average, how='left', on=['month','hour']).set_index(ef_mh.index).drop(columns=['month','hour'])

# fill with the month hourly average
for ba in list(ef_mh.columns):
    ef_mh[ba] = ef_mh[ba].fillna(mh_hold[ba])

variance_mh = pd.DataFrame()
variance_mh['stdev'] = ef_mh.std(ddof=0)

print(f'Mean Absolute percentage difference: {round(abs(((variance_mh - variance_missing) / variance_missing)).mean().values[0] *100, 2)}%')
print(f'Mean percentage difference: {round(((variance_mh - variance_missing) / variance_missing).mean().values[0] *100, 2)}%')<jupyter_output><empty_output><jupyter_text>Visualize the interpolation for a single region<jupyter_code>ba= 'GVL'
test = hourly_ef.copy()[[ba]]
test['ffill'] = test[ba].fillna(ef_ffill[ba])
test['linear'] = test[ba].fillna(ef_linear[ba])
test['MA'] = test[ba].fillna(ef_rolling[ba])
test['MA_fixed_effect'] = test[ba].fillna(ef_hourly_MA[ba])
test['monthlyTOD'] = test[ba].fillna(ef_mh[ba])


px.line(test, y=['ffill','linear','MA','MA_fixed_effect','monthlyTOD',ba])<jupyter_output><empty_output><jupyter_text>Export the interpolated data<jupyter_code>hourly_ef.to_csv('../data/processed/emission_factors/emission_factors_raw.csv')
ef_ffill.to_csv('../data/processed/emission_factors/emission_factors_ffill.csv')
ef_linear.to_csv('../data/processed/emission_factors/emission_factors_linear.csv')
ef_rolling.to_csv('../data/processed/emission_factors/emission_factors_MA.csv')
ef_hourly_MA.to_csv('../data/processed/emission_factors/emission_factors_fixed_effects.csv')
ef_mh.to_csv('../data/processed/emission_factors/emission_factors_monthhour.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1344756790903348, 'fraction_numerical': 0.020925457991156034, 'mean_word_length': 5.86449864498645, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 49, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '460929', 'license': 'MIT', 'n_tokens_mistral': 4553, 'n_tokens_neox': 4185, 'n_words': 1131, 'path': 'code/3_impute_missing_CI_data.ipynb', 'repo_name': 'grgmiller/hourly_carbon_accounting', 'size': '12664'}