Fatigue life estimation, reliability and durability are important in acquisition, maintenance and operation of vehicle systems. Fatigue life is random because of the stochastic load, the inherent variability of material properties, and the uncertainty in the definition of the S-N curve. The commonly used fatigue life estimation methods calculate the mean (not the distribution) of fatigue life under Gaussian loads using the potentially restrictive narrow-band assumption. In this paper, a general methodology is presented to calculate the statistics of fatigue life for a linear vibratory system under stationary, non-Gaussian loads considering the effects of skewness and kurtosis. The input loads are first characterized using their first four moments (mean, standard deviation, skewness and kurtosis) and a correlation structure equivalent to a given Power Spectral Density (PSD). Then, the first four moments and the correlation structure of the response process are calculated using the impulse response function of the system. Polynomial Chaos Expansion (PCE) and Karhunen-Loeve (KL) expansion are used to develop stochastic metamodels of the input and output non-Gaussian processes. Subsequently, simulated trajectories from the response stochastic metamodel are rainflow counted to obtain realizations of the fatigue life random variable based on the Miner’s damage model. Finally, the Saddlepoint Approximations (SPA) method provides the PDF and percentiles of the fatigue life. The approach can be easily extended to non-stationary loads, and with some modifications, to even non-linear vibratory systems. It can also be used to develop an accelerated testing procedure without increasing the applied load artificially. The proposed methodology is demonstrated with a quarter car example under non-Gaussian load.