A methodology is proposed for the model order reduction of finite element approximations of passive electromagnetic devices under random input conditions. In this approach, the reduced order system matrices are represented in terms of their convergent orthogonal polynomial expansions of input random variables. The coefficients of these polynomials, which are matrices, are obtained by repeated, deterministic model order reduction of finite element models generated for specific values of the input random variables. These values are chosen efficiently in a multi-dimensional grid using a Smolyak algorithm. The stochastic reduced order model is represented in the form of an augmented system which can be used for generating the desired statistics of the specific system response. The proposed method provides for significant improvement in computational efficiency over standard Monte Carlo.