<jupyter_start><jupyter_text>Stage 2 Model Selection: Summary of ResultsThis notebook generates the summary result tables and charts> Note in the paper these are figure 3 and 4.* Figure 3: Cross-validation variation in MASE across regions by horizon* Figure 4: Cross-validation variation in Prediction Interval Coverage across regions by horizon<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

import os
from os import listdir
from os.path import isfile, join

sns.set(style="whitegrid")<jupyter_output><empty_output><jupyter_text>Current working directory and file paths<jupyter_code>cwd = os.getcwd()
cwd
if cwd[-7:] != "results":
    mypath = './results/model_selection/stage2/'
    TABLE_PATH = './paper/tables/'
    FIGURE_PATH = './paper/figures/'
else:
    mypath = './model_selection/stage2/'
    TABLE_PATH = '../paper/tables/'
    FIGURE_PATH = '../paper/figures/'

result_files = [f for f in listdir(mypath) if isfile(join(mypath, f))]
results_mean = pd.DataFrame()
results_med = pd.DataFrame()
results_mean_std = pd.DataFrame()<jupyter_output><empty_output><jupyter_text>Point Estimate Results<jupyter_code>error_measures = ['smape', 'rmse', 'mase', 'coverage_80', 'coverage_95']
for metric in error_measures:
    to_read = [filename for filename in result_files if metric in filename]
    model_names = [name[:name.index('_')] for name in to_read]
    
    for filename, model_name in zip(to_read, model_names):
        df = pd.read_csv(mypath + filename, index_col=0)

        prefix = model_name + '_' + metric
        results_mean[prefix + '_mean'] = df.mean()
        results_mean[prefix  + '_std'] = df.std()
        results_med[prefix + '_med'] = df.median()
        results_med[prefix + '_iqr'] = df.quantile(0.75) - df.quantile(0.25)
        
        results_mean_std[prefix] = results_mean[prefix + '_mean'].map('{:,.2f}'.format) \
            + ' (' + results_mean[prefix  + '_std'].map('{:,.2f}'.format) + ')'
#Variability of PI Coverage across regions in the South West of England

print('Figure 4: Cross-validation variation in Prediction Interval Coverage across regions by horizon')

fig, ax = plt.subplots(2, 2, figsize=(12,8), sharey='row')
#fig.set_size_inches(w=5, h=4) 

days = [7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84]

summary_fbp = results_mean.filter(like="coverage_95").filter(like="fbp_").filter(like='mean')
ax[0][0].boxplot(x=summary_fbp.T, labels=days);
ax[0][0].axhline(y=0.95, color='red')
ax[0][0].set_ylabel('95% PI coverage')
#ax[0][0].set_xlabel('horizon (days)')
ax[0][0].set_ylim(0.88, 0.98)
ax[0][0].set_title('Prophet')
ax[0][0].xaxis.grid(False)
#ax[0][0].legend(['Prophet'], markerscale=2.0)
#ax[0][0].annotate(s='Prophet', xy=(7, 0.96), xytext=(10, 0.975))

summary_fa = results_mean.filter(like="coverage_95").filter(like="ensemble_").filter(like='mean')
ax[0][1].boxplot(x=summary_fa.T , labels=days);
ax[0][1].axhline(y=0.95, color='red')
ax[0][1].xaxis.grid(False)
ax[0][1].set_title('Ensemble')
#ax[0][1].set_xlabel('horizon (days)')
#ax[0][1].legend(['Ensemble'])

summary_fbp2 = results_mean.filter(like="coverage_80").filter(like="fbp_").filter(like='mean')
ax[1][0].boxplot(x=summary_fbp2.T, labels=days);
ax[1][0].axhline(y=0.80, color='red')
ax[1][0].xaxis.grid(False)
ax[1][0].set_ylabel('80% PI coverage')
ax[1][0].set_xlabel('horizon (days)')
#ax[1][0].legend(['Prophet'])

summary_fa2 = results_mean.filter(like="coverage_80").filter(like="ensemble_").filter(like='mean')
ax[1][1].boxplot(x=summary_fa2.T, labels=days);
ax[1][1].axhline(y=0.80, color='red')
ax[1][1].xaxis.grid(False)
ax[1][1].set_xlabel('horizon (days)')
#ax[1][1].legend(['Ensemble'])

fig.savefig(f'{FIGURE_PATH}figure4.png', dpi=300, bbox_inches='tight')
#Variability of MASE across regions in the South West of England
print('Variability of MASE by horizon (not in paper)')

fig, ax = plt.subplots(1, 2, figsize=(12,4), sharey='row')

days = [7, 14, 21, 28, 35, 42, 49, 56, 63, 70, 77, 84]

summary_fbp = results_mean.filter(like="mase").filter(like="fbp_").filter(like='mean')
ax[0].boxplot(x=summary_fbp.T, labels=days);
#ax[0][0].axhline(y=0.95, color='red')
ax[0].set_ylabel('MASE')
ax[0].set_xlabel('horizon (days)')
#ax[0].set_ylim(0.91, 0.98)
ax[0].set_title('Prophet')
ax[0].xaxis.grid(False)
#ax[0][0].legend(['Prophet'], markerscale=2.0)
#ax[0][0].annotate(s='Prophet', xy=(7, 0.96), xytext=(10, 0.975))

summary_fa = results_mean.filter(like="mase").filter(like="ensemble").filter(like='mean')
ax[1].boxplot(x=summary_fa.T , labels=days);
#ax[1].axhline(y=0.95, color='red')
ax[1].xaxis.grid(False)
ax[1].set_title('Ensemble')
ax[1].set_xlabel('horizon (days)');
#ax[0][1].legend(['Ensemble'])

region_means = results_mean.filter(like="mase").filter(like='mean').mean().sort_index()
region_std = results_mean.filter(like="mase").filter(like='std').mean().sort_index()
region_95_mean = results_mean.filter(like="coverage_95").filter(like='mean').mean().sort_index()
region_80_mean = results_mean.filter(like="coverage_80").filter(like='mean').mean().sort_index()
region_95_std = results_mean.filter(like="coverage_95").filter(like='std').mean().sort_index()
region_80_std = results_mean.filter(like="coverage_80").filter(like='std').mean().sort_index()
comparisons = list(region_means.index)
idx = [i.replace('_mase_mean', '') for i in comparisons]
df_regions = pd.DataFrame(region_means.to_numpy(), columns=['mean'])
df_regions['std'] = region_std.to_numpy()
df_regions['mean_80'] = region_80_mean.to_numpy()
df_regions['std_80'] = region_80_std.to_numpy()
df_regions['mean_95'] = region_95_mean.to_numpy()
df_regions['std_95'] = region_95_std.to_numpy()
df_regions['MASE'] = df_regions['mean'].map('{:,.2f}'.format) \
            + ' (' + df_regions['std'].map('{:,.2f}'.format) + ')'

df_regions['Coverage 80'] = df_regions['mean_80'].map('{:,.2f}'.format) \
            + ' (' + df_regions['std_80'].map('{:,.2f}'.format) + ')'

df_regions['Coverage 95'] = df_regions['mean_95'].map('{:,.2f}'.format) \
            + ' (' + df_regions['std_95'].map('{:,.2f}'.format) + ')'

df_regions.index = idx
df_regions = df_regions.drop(['mean', 'std', 'mean_80', 'std_80', 
                              'mean_95', 'std_95'], axis=1)
df_regions



from scipy.stats import norm, t
import numpy as np
alpha = 0.05
n = 67
print('Figure 3: Cross-validation variation in MASE across regions by horizon.')
fig, ax = plt.subplots(4, 2, sharex=True, figsize=(12,9))

regions = ['BNSSG', 'Cornwall', 'Devon', 'Dorset', 'Gloucestershire', 'Somerset', 'Wiltshire', 'Trust']

index = 0
for row in range(4):
    for col in range(2):
        region = regions[index]
        means = results_mean.filter(like="mase").filter(like='mean').filter(like=region)
       
        #get a consistent ordering of series (columns)
        columns = means.columns
        columns = [s.replace('_mean', '') for s in columns]
        means.columns = columns
        means = means[sorted(list(means.columns))]
        
        #plot mean by horizon
        ax[row][col].plot(means)
        
        
        #calculate upper and lower CI half width
        hw = results_mean.filter(like="mase").filter(like='std').filter(like=region)
        hw = hw[sorted(list(hw.columns))]
        hw = hw / np.sqrt(n)
        z = np.abs(t.ppf(alpha / 2, n - 1))
        
        hw *= z
        means = results_mean.filter(like="mase").filter(like='mean').filter(like=region)
        means = means[sorted(list(means.columns))]
        
        columns = hw.columns
        columns = [s.replace('_std', '') for s in columns]
        hw.columns = sorted(columns)
        means.columns = columns
        
        upper = means + hw
        lower = means - hw
        
        #col 1 = ensemble
        col1 = upper.columns[0]
        
        #col 2 = prophet
        col2 = upper.columns[1]
        
        ax[row][col].fill_between(np.arange(12), lower[col1], upper[col1], color='blue', alpha=0.15)
        
         
        ax[row][col].fill_between(np.arange(12), lower[col2], upper[col2], color='orange', alpha=0.15)
        #ax[row][col].fill_between(np.arange(12), means[col2], lower[col2], color='orange', alpha=0.2)
       
        
        ax[row][col].set_title(region)
        #ax[row][col].legend(loc='best', labels=means.columns)
        index += 1
        
ax[3][0].set_xlabel('horizon (days)')
ax[3][1].set_xlabel('horizon (days)')

ax[0][0].set_ylabel('MASE')
ax[1][0].set_ylabel('MASE')
ax[2][0].set_ylabel('MASE')
ax[3][0].set_ylabel('MASE')
ax[0][0].set_title('Bristol')

fig.legend(['Ensemble', 'Prophet'], loc='lower center',ncol=2)

fig.savefig(f'{FIGURE_PATH}figure3.png', dpi=300, bbox_inches='tight')<jupyter_output>Figure 3: Cross-validation variation in MASE across regions by horizon.
{'fraction_non_alphanumeric': 0.19079621824809204, 'fraction_numerical': 0.03952614192960474, 'mean_word_length': 4.553447185325743, 'pattern_counts': {'":': 1, '<': 12, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '160130', 'license': 'MIT', 'n_tokens_mistral': 3453, 'n_tokens_neox': 3250, 'n_words': 572, 'path': 'results/02_summary_stage_2.ipynb', 'repo_name': 'TomMonks/swast-benchmarking', 'size': '8779'}