Choosing model from Walk-Forward CV for Time Series my question is about the ‘correct’ form of analysis following the use of multiple train/test splits on time series data. Specifically, I’m using Sklearn time series split to generate 10 windows for training an XGBoost model on a sparse time series dataset (~75 time periods / rows) to not have lookahead bias. For each split, I save the MSE on the test set. 
(method similar to this link), and then average the errors over the splits as a measure of model overall performance.
I’m uncertain about the way to choose one of these models, given the different type of CV here. Do I just use the one with the lowest test error to analyze feature importance and other ‘post-modeling’ steps?
Thanks for your help. 
 A: I think this is a good chance to use nested cross validation for model selection.  If you are optimizing over hyper parameters at the same time you are doing model selection, you risk being too optimistic about out of sample error. This part of the sklearn docs does a good job of explaining nested cross validation.
Fortunately, sklearn makes it really easy to do nested cross validation with a walk forward validation scheme.  Let's set up a little problem.
I want to predict noisy observations from this function.  I'm going to use all observations before $t=20$ as train and all after $t=20$ as test (note, I'm showing the noiseless version.  The training data is noisy).

I'll use a decision tree (because I know it will do poorly) and support vector regression.  The first thing to do is set up how I will validate my models.  I'll choose TimeSeriesSplit with 5 splits as my outer validation, and TimeSeriesSplit with 5 splits as my inner validation.  Let's make some pipelines
#Make an inner and outer validation scheme
time_splitter_outer = TimeSeriesSplit(n_splits = 5)
time_splitter_inner = TimeSeriesSplit(n_splits = 5)

svr_pipe = Pipeline([
    ('scale', StandardScaler()),
    ('reg', SVR())
])

svr_params = {'reg__gamma': [0.1, 1, 2, 10]}

gs_svr = GridSearchCV(svr_pipe, param_grid = svr_params, cv = time_splitter_inner, scoring = 'neg_mean_squared_error')


tree_pipe = Pipeline([
    ('scale', StandardScaler()),
    ('reg', DecisionTreeRegressor())
])

tree_params = {'reg__max_depth': [1,2, 3]}

gs_tree = GridSearchCV(tree_pipe, param_grid = tree_params, cv = time_splitter_inner, scoring = 'neg_mean_squared_error')


Now, for the outer walk forward validation



svr_scores = cross_val_score(gs_svr, t_train, y_train, cv = time_splitter_outer, scoring = 'neg_mean_squared_error').mean()
tree_scores = cross_val_score(gs_tree, t_train, y_train, cv = time_splitter_outer, scoring = 'neg_mean_squared_error').mean()

print(svr_scores)
print(tree_scores)


We find that (unsurprisingly) the support vector machine does better.  How much better?   
SVR: -3.57
Tree: -4.03

Based on these numbers, you would choose your model.  In this case, I would choose the SVR over the tree.  Here is what the two predictions look like on the test data

So our model had a validation score of 3.57.  When predicting on the test set, we get a MSE of 3.58. Not too bad!
But, what would our score be if we had just done the grid search and not the nested CV?  It would have been 2.16!  That is way lower than our test error!  This showcases why nested CV is a really good option for this kind of validation (and for all forms of validation where we have to optimize over hyper parameters).
Here is the code to reproduce this example:
import numpy as np
from sklearn.svm import SVR
from sklearn.tree import DecisionTreeRegressor
from sklearn.model_selection import cross_val_score, GridSearchCV, TimeSeriesSplit
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
from sklearn.metrics import mean_squared_error


np.random.seed(0)
#Generate some data
t = np.linspace(0, 25, 1001)
y = np.sin(2*np.pi*t/10)*np.exp(t/12.5)
yobs = y+np.random.normal(0, 1, size = t.size)

#Keep some for testing
is_test = t>=20
t_test = t[is_test].reshape(-1,1)
y_test = y[is_test]

t_train = t[~is_test].reshape(-1,1)
y_train = y[~is_test]

#Make an inner and outer validation scheme
time_splitter_outer = TimeSeriesSplit(n_splits = 5)
time_splitter_inner = TimeSeriesSplit(n_splits = 5)

svr_pipe = Pipeline([
    ('scale', StandardScaler()),
    ('reg', SVR())
])
svr_params = {'reg__gamma': [0.01, 0.1, 1, 2, 10]}
gs_svr = GridSearchCV(svr_pipe, param_grid = svr_params, cv = time_splitter_inner, scoring = 'neg_mean_squared_error')


tree_pipe = Pipeline([
    ('scale', StandardScaler()),
    ('reg', DecisionTreeRegressor())
])
tree_params = {'reg__max_depth': [1,2, 3]}
gs_tree = GridSearchCV(tree_pipe, param_grid = tree_params, cv = time_splitter_inner, scoring = 'neg_mean_squared_error')


svr_scores = cross_val_score(gs_svr, t_train, y_train, cv = time_splitter_outer, scoring = 'neg_mean_squared_error').mean()
tree_scores = cross_val_score(gs_tree, t_train, y_train, cv = time_splitter_outer, scoring = 'neg_mean_squared_error').mean()

print('SVR:',svr_scores.round(2))
print('Tree:',tree_scores.round(2))

gs_tree.fit(t_train, y_train)
gs_svr.fit(t_train, y_train)

fig, ax = plt.subplots(dpi = 120)

plt.plot(t_test, gs_tree.predict(t_test), label = 'tree test', color = 'C0')
plt.plot(t_test, gs_svr.predict(t_test), label = 'svr test,', color = 'red')



plt.plot(t_train, gs_tree.predict(t_train), label = 'tree train', color = 'C0', alpha = 0.25)
plt.plot(t_train, gs_svr.predict(t_train), label = 'svr train', color = 'red', alpha = 0.25)

plt.plot(t[is_test], y[is_test], color = 'k', label = 'truth test' )
plt.plot(t[~is_test], y[~is_test], color = 'k', label = 'truth train', alpha = 0.25 )
plt.legend()

#Test error with selected model
print(mean_squared_error(y_test, gs_svr.predict(t_test)))

#What if we only dif GS?

gs_svr = GridSearchCV(svr_pipe, param_grid = svr_params, cv = time_splitter_inner, scoring = 'neg_mean_squared_error')

print(gs_svr.fit(t_train, y_train).best_score_)

