In many engineering applications, missing data during system identification can hinder the performance of the identified model. In this paper, a novel two-stage nonparametric framework is proposed for missing data imputation, uncertainty quantification, and its integration in system identification with reduced computational complexity. The framework does not require functional forms for both the imputation model and the identified mathematical model. Moreover, through the construction of a single imputation model, analytical expressions of predictive distributions can be given for missing entries across all missingness patterns. Furthermore, analytical expressions of the expectation and variance of distribution are provided to impute missing values and quantify uncertainty, respectively. This uncertainty is incorporated into a single mathematical model by mitigating the influence of samples with imputations during training and testing. The framework is applied to three applications, including a simulated example and two real applications on structural health monitoring and seismic attenuation modeling. Results reveal a minimum reduction of 21% in root mean squared error values, compared to those achieved by directly removing incomplete samples.