Abstract: Accurate and reliable cold, heat and electrical load forecasting is a prerequisite for optimal scheduling and efficient operation of integrated energy system. To effectively extract the linear, nonlinear, coupling and uncertainty characteristics existing among load sequences, this paper proposes a multiple linear regression（MLR）, improved complete ensemble empirical mode decomposition with adaptive noise（ICEEMDAN）, long short-time memory（LSTM） neural network, and Monte Carlo（MC） method combined with multivariate load prediction method. First, MLR models are constructed separately for cold, heat and electrical load to mine the linear features. Then, the residual components are decomposed using the ICEEMDAN method, and LSTM models are built for each load residual component in the same frequency band after the component reconstruction to realize the learning of nonlinearity and coupling. Finally, the point prediction values are acquired by summing the MLR and LSTM results. Compared with the optimal results in the reference model, the R~2 of the method improves 0.09%,0.21%,and 0.40% for cold,heat,and electrical load, respectively. In addition, to achieve an effective quantification of load uncertainty, a combination of nonparametric kernel density estimation and MC sampling is further used to obtain the prediction interval results. After the example analysis, the prediction interval coverage probability of each load is greater than the corresponding confidence level（95%,90%,85%）, and the proposed method has high prediction accuracy and reliability.