Multiple Comparisons

This setup code is required to run in an IPython notebook

In [1]:
import warnings
warnings.simplefilter('ignore')

%matplotlib inline
import seaborn
seaborn.mpl.rcParams['figure.figsize'] = (10.0, 6.0)
seaborn.mpl.rcParams['savefig.dpi'] = 90
seaborn.mpl.rcParams['font.family'] = 'serif'
seaborn.mpl.rcParams['font.size'] = 14

# Reproducability
import numpy as np
np.random.seed(23456)
# Common seed used throughout
seed = np.random.randint(0,2**31-1)

The multiple comparison procedures all allow for examining aspects of superior predictive ability. There are three available:

  • SPA - The test of Superior Predictive Ability, also known as the Reality Check (and accessible as RealityCheck) or the bootstrap data snooper, examines whether any model in a set of models can outperform a benchmark.
  • StepM - The stepwise multiple testing procedure uses sequential testing to determine which models are superior to a benchmark.
  • MCS - The model confidence set which computes the set of models which with performance indistinguishable from others in the set.

All procedures take losses as inputs. That is, smaller values are preferred to larger values. This is common when evaluating forecasting models where the loss function is usually defined as a positive function of the forecast error that is increasing in the absolute error. Leading examples are Mean Square Error (MSE) and Mean Absolute Deviation (MAD).

The test of Superior Predictive Ability (SPA)

This procedure requires a $t$-element array of benchmark losses and a $t$ by $k$-element array of model losses. The null hypothesis is that no model is better than the benchmark, or

$$ H_0: \max_i E[L_i] \geq E[L_{bm}] $$

where $L_i$ is the loss from model $i$ and $L_{bm}$ is the loss from the benchmark model.

This procedure is normally used when there are many competing forecasting models such as in the study of technical trading rules. The example below will make use of a set of models which are all equivalently good to a benchmark model and will serve as a size study.

Study Design

The study will make use of a measurement error in predictors to produce a large set of correlated variables that all have equal expected MSE. The benchmark will have identical measurement error and so all models have the same expected loss, although will have different forecasts.

The first block computed the series to be forecast.

In [2]:
from numpy.random import randn
import statsmodels.api as sm
t = 1000
factors = randn(t,3)
beta = np.array([1,0.5,0.1])
e = randn(t)
y = factors.dot(beta)

The next block computes the benchmark factors and the model factors by contaminating the original factors with noise. The models are estimated on the first 500 observations and predictions are made for the second 500. Finally, losses are constructed from these predictions.

In [3]:
# Measurement noise
bm_factors = factors + randn(t,3)
# Fit using first half, predict second half
bm_beta = sm.OLS(y[:500],bm_factors[:500]).fit().params
# MSE loss
bm_losses = (y[500:] - bm_factors[500:].dot(bm_beta))**2.0
# Number of models
k = 500
model_factors = np.zeros((k,t,3))
model_losses = np.zeros((500,k))
for i in range(k):
    # Add measurement noise
    model_factors[i] = factors + randn(1000,3)
    # Compute regression parameters
    model_beta  = sm.OLS(y[:500],model_factors[i,:500]).fit().params
    # Prediction and losses
    model_losses[:,i] = (y[500:] - model_factors[i,500:].dot(model_beta))**2.0

Finally the SPA can be used. The SPA requires the losses from the benchmark and the models as inputs. Other inputs allow the bootstrap sued to be changed or for various options regarding studentization of the losses. compute does the real work, and then pvalues contains the probability that the null is true given the realizations.

In this case, one would not reject. The three p-values correspond to different re-centerings of the losses. In general, the consistent p-value should be used. It should always be the case that

$$ lower \leq consistent \leq upper .$$

See the original papers for more details.

In [4]:
from arch.bootstrap import SPA
spa = SPA(bm_losses, model_losses)
spa.seed(seed)
spa.compute()
spa.pvalues
Out[4]:
lower         0.520
consistent    0.723
upper         0.733
dtype: float64

The same blocks can be repeated to perform a simulation study. Here I only use 100 replications since this should complete in a reasonable amount of time. Also I set reps=250 to limit the number of bootstrap replications in each application of the SPA (the default is a more reasonable 1000).

In [5]:
# Save the pvalues
pvalues = []
b = 100
seeds = np.random.randint(0, 2**31 - 1, b)
# Repeat 100 times
for j in range(b):
    if j % 10 == 0:
        print(j)
    factors = randn(t,3)
    beta = np.array([1,0.5,0.1])
    e = randn(t)
    y = factors.dot(beta)


    # Measurement noise
    bm_factors = factors + randn(t,3)
    # Fit using first half, predict second half
    bm_beta = sm.OLS(y[:500],bm_factors[:500]).fit().params
    # MSE loss
    bm_losses = (y[500:] - bm_factors[500:].dot(bm_beta))**2.0
    # Number of models
    k = 500
    model_factors = np.zeros((k,t,3))
    model_losses = np.zeros((500,k))
    for i in range(k):
        model_factors[i] = factors + randn(1000,3)
        model_beta  = sm.OLS(y[:500],model_factors[i,:500]).fit().params
        # MSE loss
        model_losses[:,i] = (y[500:] - model_factors[i,500:].dot(model_beta))**2.0
    # Lower the bootstrap replications to 250
    spa = SPA(bm_losses, model_losses, reps = 250)
    spa.seed(seeds[j])
    spa.compute()
    pvalues.append(spa.pvalues)
0
10
20
30
40
50
60
70
80
90

Finally the pvalues can be plotted. Ideally they should form a $45^o$ line indicating the size is correct. Both the consistent and upper perform well. The lower has too many small p-values.

In [6]:
import pandas as pd

pvalues = pd.DataFrame(pvalues)
for col in pvalues:
    values = pvalues[col].values
    values.sort()
    pvalues[col] = values
# Change the index so that the x-values are between 0 and 1
pvalues.index = np.linspace(0.005,.995,100)
fig = pvalues.plot()

Power

The SPA also has power to reject then the null is violated. The simulation will be modified so that the amount of measurement error differs across models, and so that some models are actually better than the benchmark. The p-values should be small indicating rejection of the null.

In [7]:
# Number of models
k = 500
model_factors = np.zeros((k,t,3))
model_losses = np.zeros((500,k))
for i in range(k):
    scale = ((2500.0 - i) / 2500.0)
    model_factors[i] = factors + scale * randn(1000,3)
    model_beta  = sm.OLS(y[:500],model_factors[i,:500]).fit().params
    # MSE loss
    model_losses[:,i] = (y[500:] - model_factors[i,500:].dot(model_beta))**2.0

spa = SPA(bm_losses, model_losses)
spa.seed(seed)
spa.compute()
spa.pvalues
Out[7]:
lower         0.0
consistent    0.0
upper         0.0
dtype: float64

Here the average losses are plotted. The higher index models are clearly better than the lower index models -- and the benchmark model (which is identical to model.0).

In [8]:
model_losses = pd.DataFrame(model_losses,columns=['model.' + str(i) for i in range(k)])
avg_model_losses = pd.DataFrame(model_losses.mean(0), columns=['Average loss'])
fig = avg_model_losses.plot(style=['o'])

Stepwise Multiple Testing (StepM)

Stepwise Multiple Testing is similar to the SPA and has the same null. The primary difference is that it identifies the set of models which are better than the benchmark, rather than just asking the basic question if any model is better.

In [9]:
from arch.bootstrap import StepM
stepm = StepM(bm_losses, model_losses)
stepm.compute()
print('Model indices:')
print([model.split('.')[1] for model in stepm.superior_models])
Model indices:
['106', '115', '117', '152', '156', '157', '158', '169', '186', '187', '192', '197', '214', '215', '219', '228', '235', '237', '238', '244', '246', '248', '252', '253', '254', '257', '261', '262', '263', '266', '272', '275', '279', '280', '281', '282', '286', '291', '294', '298', '299', '300', '305', '306', '310', '312', '315', '316', '318', '325', '326', '328', '329', '330', '332', '333', '335', '336', '338', '340', '341', '342', '344', '348', '349', '350', '351', '352', '353', '354', '356', '357', '359', '360', '362', '363', '364', '365', '366', '367', '368', '370', '371', '372', '373', '374', '375', '377', '378', '379', '380', '382', '383', '385', '386', '387', '388', '389', '390', '391', '392', '393', '394', '395', '398', '399', '400', '401', '402', '403', '404', '405', '406', '407', '408', '410', '411', '412', '413', '414', '417', '419', '420', '421', '422', '423', '424', '425', '426', '427', '428', '429', '431', '432', '433', '434', '435', '436', '437', '438', '439', '440', '441', '442', '443', '444', '445', '447', '448', '449', '450', '451', '453', '454', '455', '456', '457', '458', '459', '460', '461', '462', '463', '464', '465', '466', '467', '468', '469', '470', '471', '473', '474', '475', '476', '477', '478', '479', '480', '481', '482', '483', '484', '485', '486', '487', '488', '489', '490', '491', '492', '493', '494', '495', '496', '497', '498', '499', '73']
In [10]:
better_models = pd.concat([model_losses.mean(0),model_losses.mean(0)],1)
better_models.columns = ['Same or worse','Better']
better = better_models.index.isin(stepm.superior_models)
worse = np.logical_not(better)
better_models.loc[better,'Same or worse'] = np.nan
better_models.loc[worse,'Better'] = np.nan
fig = better_models.plot(style=['o','s'], rot=270)

The Model Confidence Set

The model confidence set takes a set of losses as its input and finds the set which are not statistically different from each other while controlling the familywise error rate. The primary output is a set of p-values, where models with a pvalue above the size are in the MCS. Small p-values indicate that the model is easily rejected from the set that includes the best.

In [11]:
from arch.bootstrap import MCS
# Limit the size of the set
losses = model_losses.iloc[:,::20]
mcs = MCS(losses, size=0.10)
mcs.compute()
print('MCS P-values')
print(mcs.pvalues)
print('Included')
included = mcs.included
print([model.split('.')[1] for model in included])
print('Excluded')
excluded = mcs.excluded
print([model.split('.')[1] for model in excluded])
MCS P-values
            Pvalue
Model name        
model.60     0.000
model.80     0.001
model.40     0.002
model.140    0.002
model.20     0.005
model.100    0.005
model.120    0.013
model.0      0.013
model.220    0.025
model.160    0.128
model.240    0.130
model.260    0.146
model.200    0.146
model.180    0.403
model.320    0.460
model.420    0.483
model.400    0.717
model.360    0.866
model.340    0.866
model.280    0.866
model.460    0.866
model.380    0.866
model.300    0.866
model.480    0.866
model.440    1.000
Included
['160', '180', '200', '240', '260', '280', '300', '320', '340', '360', '380', '400', '420', '440', '460', '480']
Excluded
['0', '100', '120', '140', '20', '220', '40', '60', '80']
In [12]:
status = pd.DataFrame([losses.mean(0),losses.mean(0)], index=['Excluded','Included']).T
status.loc[status.index.isin(included), 'Excluded'] = np.nan
status.loc[status.index.isin(excluded), 'Included'] = np.nan
fig = status.plot(style=['o','s'])