A novel algorithm is proposed for simulating univariate non-Gaussian nonstationary processes (NNP) with the specified evolutionary power spectral density (EPSD)/nonstationary auto-correlation function (NACF) and first four-order time-varying marginal moments (TVMMs). The sample realizations of the target NNP are generated as the outputs from a specific time-varying auto-regressive (TVAR) model via filtering the non-Gaussian and nonstationary white noise inputs. These white noise inputs are also non-Gaussian and nonstationary, and their first four-order TVMMs are predetermined using an approach developed herein according to the specified EPSD/NACF and first four-order TVMMs of the outputs. The conventional Johnson transformation is updated to accommodate the nonstationary cases for producing desired white noise inputs. This algorithm is developed from the linear filtering method (LFM), and inherits the simplicity and high efficiency from LFM. It fills the gaps in LFM-based algorithms for simulating NNP. Two numerical examples, i.e., a ground motion acceleration and a downburst velocity, are presented to fully demonstrate the capabilities of the proposed algorithm by comparing the simulation statistics with the targets.