Roman Kh is correct to warn you against ever using stepwise approaches. One of the best discussions of their pitfalls is Peter Flom's paper Stop Using Stepwise
http://www.lexjansen.com/pnwsug/2008/DavidCassell-StoppingStepwise.pdf
That said, every statistician and their brother has a paper or approach to variable selection – they are legion. To your point, these are all focused on selection with a single response function. I am not aware of anyone who has developed algorithms specifically for use with multiple dependent variables and would be happy to be told that this is incorrect, someone, somewhere has a protocol.
Methodological solutions follow demand and, if there is no demand, then no one will bother. To date, there would appear to be, at best, limited demand for selection routines with multiple response functions. I really do not understand why nearly all modeling projects insist on choosing a single response function when multiple functions would give a much better, more informative and insightful answer. There are many possible reasons for this but, in my view, the leading explanations would have to include a deeply engrained bias in favor of “Occam’s Razor-like,” single response models; the paucity of training in the use and interpretation of multiple response models as well as the inevitable consequences of our cognitive limitations in “bounded rationality.” This is true despite the fact that all of the major statistical packages offer MANOVA or canonical correlation routines. What they all lack is a “LASSO-like” algorithm for multiple DVs and large numbers of candidate features.
An informative exception to these observations is a paper by Grice and Iwasaki which compares ANOVA with MANOVA, discussing the advantages and pitfalls of each in the context of hypothesis-testing, inference and interpretation. Note that they do not address your specific issue concerned with variable selection.
http://psychology.okstate.edu/faculty/jgrice/psyc6813/Grice_Iwasaki_AMR.pdf
This paper raises a fundamental issue which the OP hasn’t addressed: the objectives of the model. Is it to be used for hypothesis-testing and inference or black box prediction as in a machine learning problem? These really are independent challenges with differing solutions in large part as a function of the amount of information under analysis. For relatively small amounts of data, classic inferential methods are realistic. If one is faced with large amounts of information containing many, even massive quantities of candidate predictors, then the classic approaches break down.
Given this, what are the limiting cases for variable selection with multiple dependent variables? Of course, one always has the option of combining the multiple DVs into an a priori composite. In this instance, the variable selection process would be the same as for any single response function.
When modeling truly multiple DVs the simplest and most obvious example would be to have such a finite amount of information possessing so few possible features that variable selection becomes moot, permitting the ready fitting of a canonical correlation or MANOVA as in the Grice and Iwasaki paper. This case would be consistent for use with a PhD dissertation or paper employing careful, classic hypothesis-tests.
For the more likely case where there are a large number of candidate predictors – making a variable selection step unavoidable -- a brute force solution might be to fit a separate selection process for each dependent variable. This approach should not be recommended and is flawed in that it ignores the linear combinations or composites that are inherent to a truly multivariate approach and begs the question of how a rigorous and final variable selection process would work.
It would appear to be the case that classic multivariate statistics and analysis does not offer a solution to the problem of variable selection for multiple response functions with large, even massive, numbers of candidate predictors and/or “big” data. In my view, this necessitates employing approximating workarounds that involve extensions of Breiman’s random forests routine. Breiman discussed using RFs as a variable selection method but never said that you could not employ a multivariate tool other than CART as the engine driving the algorithm. Breiman’s classic approach to RFs was limited in that it was developed in the 90s for only a few thousand candidate predictors on a single CPU. In the applied world of today, access to massively parallel platforms (MPP) for crunching massive amounts of data as well as “divide and conquer” routines means that one is no longer limited to his classic solution. For a discussion of “D&C” routines, see this paper by Chen and Minge A Split-and-Conquer Approach for Analysis of Extraordinarily Large Data
http://dimacs.rutgers.edu/TechnicalReports/TechReports/2012/2012-01.pdf
One example of how a “D&C” approach might work on an MPP for multiple response functions and thousands or hundreds of thousands (or more) candidate predictors – a common challenge with unstructured information -- would be to plug in MANOVA or canonical correlation, run millions of “mini-models” and aggregate the output on the back end to obtain both ensemble predictions as well as a ranking of truly multivariate variable relative importance. This could be done in a few hours on an MPP of reasonable size. Given the approximating nature of this approach, the modeler is forced to give up any notion of finding a final, reduced or fixed set of mathematically unique predictors. Note, however, that this would facilitate the elimination of large numbers of candidate variables.
At this point, the question becomes one of whether or not this solution is an end in itself -- is the objective prediction or inference? If it is prediction, this this could be the end product and retaining the results from the millions of mini-models would enable their later use in scoring new data. If inference is the objective, then it’s not the end of the analysis and further refinement of the variables in additional stages of modeling would further reduce the variables as well as eliminate the inevitable redundancies and pure linear combinations hidden in the rankings. At this stage of the development of D&C routines, there don’t seem to be any good answers as to how best to pursue additional stages of inferential modeling.
Anyway, these are just a few thoughts. Hope they’re helpful.