<jupyter_start><jupyter_text>UtilsFunctions to add date components as columns, converts a series to a supervised machine learning problem, split train and test, rmse, smape and printing scores of models<jupyter_code>import math 
import re


# From fast.ai 0.7
def add_datepart(df, fldname, drop=True, time=False):
    "Helper function that adds columns relevant to a date."
    df = df.copy()
    fld = df[fldname]
    fld_dtype = fld.dtype
    if isinstance(fld_dtype, pd.core.dtypes.dtypes.DatetimeTZDtype):
        fld_dtype = np.datetime64

    if not np.issubdtype(fld_dtype, np.datetime64):
        df[fldname] = fld = pd.to_datetime(fld, infer_datetime_format=True)
    targ_pre = re.sub('[Dd]ate$', '', fldname)
    attr = ['Year', 'Month', 'Week', 'Day', 'Dayofweek', 'Dayofyear',
            'Is_month_end', 'Is_month_start', 'Is_quarter_end', 'Is_quarter_start', 'Is_year_end', 'Is_year_start']
    if time: attr = attr + ['Hour', 'Minute', 'Second']
    for n in attr: df[targ_pre + n] = getattr(fld.dt, n.lower())
    df[targ_pre + 'Elapsed'] = fld.astype(np.int64) // 10 ** 9
    if drop: df.drop(fldname, axis=1, inplace=True)
    return df
        
        
def timeseries_to_supervised(df, lag=1, dropnan=True):
    
    df = pd.DataFrame(df)
    columns = [df.shift(i).rename(columns=lambda x: f"{x}(-{i})") for i in range(1, lag+1)]
    columns.insert(0, df)
    df = pd.concat(columns, axis=1)
    
    if dropnan:
        df.dropna(inplace=True)
    else:
        df = df[lag:] #skip nan at the start
        df.fillna(-1, inplace=True)
        
    return df


def split_vals(a,n): return a[:n].copy(), a[n:].copy()

def rmse(x,y): return math.sqrt(((x-y)**2).mean())

def print_score(m):
    res = [rmse(m.predict(X_train), y_train), rmse(m.predict(X_valid), y_valid),
                m.score(X_train, y_train), m.score(X_valid, y_valid)]
    if hasattr(m, 'oob_score_'): res.append(m.oob_score_)
    print(res)

def smape(predicted, actual):
    dividend= np.abs(np.array(actual) - np.array(predicted))
    denominator = np.array(actual) + np.array(predicted)
    
    return 2 * np.mean(np.divide(dividend, denominator, out=np.zeros_like(dividend), where=denominator!=0, casting='unsafe'))<jupyter_output><empty_output><jupyter_text>Load merged air quality and<jupyter_code>## Load 

donggaocun_df = pd.read_feather(DATA_PATH + 'interim/aq_station_merged/donggaocun_aq_merged_weather.feather')
donggaocun_df.iloc[9996:, :]
donggaocun_df.head()<jupyter_output><empty_output><jupyter_text>Checking missing data<jupyter_code>aq_missing_data = donggaocun_df[['PM2.5', 'PM10', 'NO2', 'CO', 'O3', 'SO2']].copy()
aq_missing_data = aq_missing_data.rename(columns={'NO2': 'Nitrogen Dioxide', 'CO': 'Carbon Monoxide', 'O3': 'Ozone', 'SO2': 'Sulfur Dioxide'})
msno.matrix(aq_missing_data)
huairou = pd.read_feather(DATA_PATH + 'interim/aq_station_merged/huairou_aq_merged_weather.feather')
fig, ax = plt.subplots(figsize=(9, 7))
donggaocun_df.set_index('utc_datetime')['PM2.5'].loc['March 31, 2017'].plot(ax=ax, label='Donggaocun')
other_aq.set_index('utc_datetime')['PM2.5'].loc['March 31, 2017'].plot(ax=ax, label='Huairou')
ax.set_title('PM 2.5 fluctuations on March 31, 2017',  fontsize=18)
ax.set_xlabel('Hours',  fontsize=18)
ax.set_ylabel('PM 2.5 microgram/cubic meter',  fontsize=18)
ax.legend(loc='best')
fig.savefig('PM25_day.png')<jupyter_output><empty_output><jupyter_text>Imputer<jupyter_code>from fancyimpute import KNN
from fancyimpute import IterativeImputer


def impute_values(df, imputer=None, columns_to_impute=None):
    
    if imputer and (columns_to_impute is not None):
            df_remaining = df.drop(columns_to_impute, axis=1)
            df_to_impute = df[columns_to_impute]
            imputed_values = imputer.fit_transform(df_to_impute)
            df_imputed = pd.DataFrame(imputed_values, columns=columns_to_impute, index=df_remaining.index)
            df = pd.concat([df_remaining, df_imputed], axis=1)
    return df
    
columns_to_impute = donggaocun_df.columns[donggaocun_df.columns.str.contains('meo|grid')]
columns_to_impute = columns_to_impute[~columns_to_impute.str.contains('weather')]
imputed_df = impute_values(donggaocun_df, KNN(5), columns_to_impute)<jupyter_output>Imputing row 1/10094 with 0 missing, elapsed time: 17.928
Imputing row 101/10094 with 0 missing, elapsed time: 17.928
Imputing row 201/10094 with 0 missing, elapsed time: 17.930
Imputing row 301/10094 with 0 missing, elapsed time: 17.931
Imputing row 401/10094 with 0 missing, elapsed time: 17.932
Imputing row 501/10094 with 0 missing, elapsed time: 17.932
Imputing row 601/10094 with 0 missing, elapsed time: 17.933
Imputing row 701/10094 with 0 missing, elapsed time: 17.934
Imputing row 801/10094 with 0 missing, elapsed time: 17.934
Imputing row 901/10094 with 15 missing, elapsed time: 17.935
Imputing row 1001/10094 with 15 missing, elapsed time: 17.941
Imputing row 1101/10094 with 0 missing, elapsed time: 17.945
Imputing row 1201/10094 with 0 missing, elapsed time: 17.945
Imputing row 1301/10094 with 0 missing, elapsed time: 17.946
Imputing row 1401/10094 with 0 missing, elapsed time: 17.947
Imputing row 1501/10094 with 0 missing, elapsed time: 17.947
Imputing row 1601/10094 with 0 mis[...]<jupyter_text>Add datetime components and dummy weather features<jupyter_code>processed_df = pd.get_dummies(imputed_df)
add_datepart(processed_df, fldname='utc_datetime',drop=False, time=True)
pd.read_feather(DATA_PATH + 'processed/donggaocun_processed.feather').columns<jupyter_output><empty_output><jupyter_text>Generate X_train, y_train, X_valid and y_valid<jupyter_code>def generate_train_test_data(df, target_column, split_ratio=None, split_date=None, lag_observations=1, dropnan=True):
    
    supervised_df = timeseries_to_supervised(df, lag_observations, dropnan=dropnan)

    if split_ratio is not None:
        split_index = int(len(supervised_df) * split_ratio)
    elif split_date is not None:
        date_index = supervised_df.index.get_loc(split_date)
        if type(date_index) == slice:
            split_index = date_index.start
        else:
            split_index = date_index
    else:
        raise("Must specify split ratio or split date")
    
    non_lagged_columns = list(supervised_df.columns[~supervised_df.columns.str.contains('\(.*\)')])
    non_lagged_columns.append(target_column)
    X = supervised_df.drop(non_lagged_columns, axis=1)
    y = supervised_df[target_column]
    
    
    X_train, X_valid = split_vals(X, split_index)
    y_train, y_valid = split_vals(y, split_index)
       
    
    return X_train, y_train, X_valid, y_valid
    
    
processed_df.set_index('utc_datetime', inplace=True)
processed_df.reset_index().to_feather(DATA_PATH + f"processed/donggaocun_processed.feather")
target_column='O3'
X_train, y_train, X_valid, y_valid = generate_train_test_data(processed_df, split_ratio=0.8, target_column=target_column, lag_observations=15, dropnan=True)<jupyter_output><empty_output><jupyter_text>Test initial model<jupyter_code>rf = RandomForestRegressor(n_estimators=100, n_jobs=-1, oob_score=True)
rf.fit(X_train, y_train)
print_score(rf)<jupyter_output><empty_output><jupyter_text>Test optimal number of lagging hours per air pollutant<jupyter_code>def find_optimal_lagging_hours(df, target_column, split_ratio=0.8):
    rmse_results = []
    smape_results = []
    for hours_lagged in range(1, 24):

        X_train, y_train, X_valid, y_valid = generate_train_test_data(df, split_ratio=split_ratio, target_column=target_column, lag_observations=hours_lagged, dropnan=True)        
        
        rf = RandomForestRegressor(n_estimators=100, n_jobs=-1)
        rf.fit(X_train, y_train)
        predictions = rf.predict(X_valid)
        rmse_results.append(rmse(predictions, y_valid))
        smape_results.append(smape(predictions, y_valid))
    
    return rmse_results, smape_results
    
pm25_rmse_results, pm25_smape_results = find_optimal_lagging_hours(processed_df, 'PM2.5')
pd.Series(pm25_rmse_results).sort_values().head()
pd.Series(pm25_smape_results).sort_values().head()
pd.Series(pm25_rmse_results).plot()
pd.Series(pm25_smape_results).plot()<jupyter_output><empty_output><jupyter_text>**Select n=16 as number of lagged observations for PM2.5**<jupyter_code>pm10_rmse_results, pm10_smape_results = find_optimal_lagging_hours(processed_df, 'PM10')
pd.Series(pm10_rmse_results).sort_values().head()
pd.Series(pm10_rmse_results).sort_values().plot()
pd.Series(pm10_smape_results).sort_values().min(), pd.Series(pm10_smape_results).sort_values().min()<jupyter_output><empty_output><jupyter_text>**Select n=13 as number of lagged observations for PM10. PM10 SMAPE has low variance for lagged observations**<jupyter_code>o3_rmse_results, o3_smape_results = find_optimal_lagging_hours(processed_df, 'O3')
pd.Series(o3_rmse_results).sort_values().head()
pd.Series(o3_smape_results).sort_values().head(15)<jupyter_output><empty_output><jupyter_text>**Select n=10 as number of lagged observations for O3.**  Loading and visualizing XGBoost and LightGBM results<jupyter_code>pm25_actual = pd.read_csv(DATA_PATH + 'processed/pm25_actual.csv', header=None, index_col=0, squeeze=True)
pm10_actual = pd.read_csv(DATA_PATH + 'processed/pm10_actual.csv', header=None, index_col=0, squeeze=True)
O3_actual = pd.read_csv(DATA_PATH + 'processed/O3_actual.csv', header=None, index_col=0, squeeze=True)<jupyter_output><empty_output><jupyter_text>**XGBoost predictions**<jupyter_code>pm25_xgb_predictions = np.load(DATA_PATH + 'processed/pm25_xg_predictions.npy')
pm25_xgb_predictions = pd.Series(pm25_xgb_predictions, index=pm25_actual.index)

pm10_xgb_predictions = np.load(DATA_PATH + 'processed/pm10_xg_predictions.npy')
pm10_xgb_predictions = pd.Series(pm10_xgb_predictions, index=pm10_actual.index)

O3_xgb_predictions = np.load(DATA_PATH + 'processed/O3_xg_predictions.npy')
O3_xgb_predictions = pd.Series(O3_xgb_predictions, index=O3_actual.index)
# pm25_actual.index.min(), pm25_actual.index.max()
# pm25_actual.shape
# pd.date_range(pm25_actual.index.min(), pm25_actual.index.max(), freq='H')<jupyter_output><empty_output><jupyter_text>**Plotting XGBoost predictions vs. Actual values**<jupyter_code>fig, ax = plt.subplots(figsize=(15, 5))

index = pd.to_datetime(pm25_actual.index)

xaxis_formatter = DateFormatter("%m %d %Y")
ax.xaxis.set_major_formatter(xaxis_formatter)


ax.plot(index, pm25_actual, label='Actual', color='blue')
ax.plot(index, pm25_xgb_predictions, label='Predicted', color='red')

ax.set_title('PM 2.5 XGBoost Predictions (3 months)')
ax.legend()

fig.autofmt_xdate()

fig.savefig('pm25_xgb.png')


fig, ax = plt.subplots(figsize=(15, 5))

index = pd.to_datetime(pm10_actual.index)

xaxis_formatter = DateFormatter("%m %d %Y")
ax.xaxis.set_major_formatter(xaxis_formatter)


ax.plot(index, pm10_actual, label='Actual', color='blue')
ax.plot(index, pm10_xgb_predictions, label='Predicted', color='red')

ax.set_title('PM 10 XGBoost Predictions (3 months)')
ax.legend()

fig.autofmt_xdate()

fig.savefig('pm10_xgb.png')

fig, ax = plt.subplots(figsize=(15, 5))

plt.plot(O3_actual, label='O3 Actual')
plt.plot(O3_xgb_predictions, label='Predictions')
ax.set_title('O3 XGBoost Predictions (2 months)')
fig.savefig('O3_xgb.png')<jupyter_output><empty_output><jupyter_text>**LightGBM predictions**<jupyter_code>pm25_lgb_predictions = np.load(DATA_PATH + 'processed/pm25_lgb_predictions.npy')
pm10_lgb_predictions = np.load(DATA_PATH + 'processed/pm10_lgb_predictions.npy')
O3_lgb_predictions = np.load(DATA_PATH + 'processed/O3_lgb_predictions.npy')
pd.to_datetime(pm25_actual.index)
from matplotlib.dates import DateFormatter

fig, ax = plt.subplots(figsize=(15, 5))

index = pd.to_datetime(pm25_actual.index)

xaxis_formatter = DateFormatter("%m %d %H")
ax.xaxis.set_major_formatter(xaxis_formatter)


ax.plot(index, pm25_actual, label='Actual', color='blue')
ax.plot(index, pm25_lgb_predictions, label='Predicted', color='red')

ax.set_title('PM 2.5 LightGBM Predictions (3 months)')
ax.legend()

fig.autofmt_xdate()


fig.savefig('pm25_lgb.png')
pm10_actual
fig, ax = plt.subplots(figsize=(15, 5))

index = pd.to_datetime(pm10_actual.index)

xaxis_formatter = DateFormatter("%m %d %H")
ax.xaxis.set_major_formatter(xaxis_formatter)

ax.plot(index, pm10_actual, label='Actual', color='blue')
ax.plot(index, pm10_lgb_predictions, label='Predicted', color='red')

ax.set_title('PM 10 LightGBM Predictions (2 months)')
ax.legend()

fig.autofmt_xdate()

fig.savefig('pm10_lgb.png')

fig, ax = plt.subplots(figsize=(20, 5))

plt.plot(O3_actual)
plt.plot(O3_lgb_predictions)
ax.set_title('O3 LightGBM Predictions (2 months)')
fig.savefig('O3_lgb.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1281222913875975, 'fraction_numerical': 0.040737530533448904, 'mean_word_length': 5.9015769439913, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '595943', 'license': 'MIT', 'n_tokens_mistral': 5224, 'n_tokens_neox': 4807, 'n_words': 840, 'path': 'notebooks/Decision_Tree_Ensembles.ipynb', 'repo_name': 'mjjimenez/kdd_air_quality_forecasting', 'size': '12691'}