We estimate the latent factors driving the non-normal dependence in high dimensional panel data using Higher-order multi-cumulant Factor Analysis (HFA). This new approach consists of an eigenvalue ratio test to select the number of non-Gaussian factors, and uses alternating regressions to estimate both the Gaussian and non-Gaussian factors. Simulation results confirm that the HFA estimators improve the accuracy of factor selection and factor estimation as compared to approaches using principal or independent component analysis. The empirical usefulness of the HFA approach is shown in an application to forecasting the S&P 500 equity premium.