The Shapely value originates from the cooperative game theory and has been widely used in solving machine learning problems due to its high interpretability and consistency. The typical Shapley value evaluates the importance score of a feature as its average marginal contribution to a fully parameterized model under all possible feature combinations; as a result, its value includes the influences from both selected features and unselected ones. To better separate the corresponding sources, it is suggested to decompose the Shapley value into high-order interaction effect components such that the influences from each individual feature can be effectively untangled. A feature's contribution is decomposed into distinct interaction components, and each component corresponds to a joint contribution resulting from a particular feature combination. The feature ranking is therefore evaluated with respect to the selected feature subset and calculated as the total incremental contribution by summing up its corresponding decomposed interaction values. In this study, a computationally efficient model-dependent greedy search algorithm on the high-order interaction components is proposed to solve an optimal subset selection problem with hundreds of time-lagged interrelated input features. Our algorithm extends a recently developed low-order polynomial time method for calculating the interaction component values. The empirical analysis demonstrates that the proposed method always outperforms other methods that are based on the typical Shapley value, the gain or the split count criteria, in terms of in-sample representativeness and out-of-sample forecasting performance for handling a problem with hundreds of time-lagged input features.