Focusing on identification, support recovery, and parameter estimation of multivariate ARMA systems with exogenous inputs (ARMAX) is considered. It is assumed that sparsity appears not only in the unknown parameter matrix of the system inputs and outputs but also in the coefficients of system noises, and with exogenous inputs the system admits feedback control and thus the observation of the ARMAX system may be nonstationary. A two-stage sparse identification algorithm is introduced. First, the extended least squares (ELS) algorithm is applied to obtain an initial estimate for the unknown parameter matrix of the ARMAX system. Second, a convex optimization criterion is introduced based on norm with regularization, where adaptive weights generated from ELS algorithms are adopted in the optimization variables in the term. Third, it is proved that by optimizing the criterion, the sets of the zero and the nonzero entries of the parameter matrix, i.e., the support set of the ARMAX system, can be correctly identified with a finite number of observations, and estimates of the nonzero parameters converge to the true values with probability one. Fourth, the strictly positive realness (SPR) condition required by the ELS algorithm is relaxed by an overparametrization technique. The performance of the algorithm is testified through simulation examples.