Treatment effect estimates in regression discontinuity (RD) designs are often sensitive to the choice of bandwidth and polynomial order, the two important ingredients of widely used local regression methods. While Imbens and Kalyanaraman (2012) and Calonico, Cattaneo and Titiunik (2014) provide guidance on bandwidth, the sensitivity to polynomial order still poses a conundrum to RD practitioners. It is understood in the econometric literature that applying the argument of bias reduction does not help resolve this conundrum, since it would always lead to preferring higher orders. We therefore extend the frameworks of Imbens and Kalyanaraman (2012) and Calonico,Cattaneo and Titiunik (2014) and use the asymptotic mean squared error of the local regression RD estimator as the criterion to guide polynomial order selection. We show in Monte Carlo simulations that the proposed order selection procedure performs well, particularly in large sample sizes typically found in empirical RD applications. This procedure extends easily to fuzzy regression discontinuity and regression kink designs.