Abstract Industrial load takes a big portion of the total electricity demand. Skilled probabilistic industrial load forecasts allow for optimally exploiting energy resources, managing the reserves, and market bidding, which are beneficial to transmission and distribution system operators and their industrial customers. Despite its importance, industrial load forecasting has never been a popular subject in the literature. Most existing methods operate on the active power alone, partially or totally neglecting the reactive power. This paper proposes a multivariate approach to probabilistic industrial load forecasting, which addresses active and reactive power simultaneously. The proposed method is based on a two-level procedure, which consists of generating probabilistic forecasts individually for active and reactive power through univariate probabilistic models, and combining these forecasts in a multivariate approach based on a multivariate quantile regression model. The procedure to estimate the parameters of the multivariate quantile regression model is posed in this paper under a linear programming problem, to facilitate the convergence to the optimal solution. The proposed method is validated using actual load data collected at an Italian factory, under comparison with several probabilistic benchmarks. The proposed multivariate method enhances the skill of forecasts by 6% to 13.5%, with respect to univariate benchmarks.