Bias correction methods are used to calibrate climate model outputs with respect to observational records. The goal is to ensure that statistical features (such as means and variances) of climate simulations are coherent with observations. In this article, a multivariate stochastic bias correction method is developed based on optimal transport. Bias correction methods are usually defined as transfer functions between random variables. We show that such transfer functions induce a joint probability distribution between the biased random variable and its correction. The optimal transport theory allows us constructing a joint distribution that minimizes an energy spent in bias correction. This extends the classical univariate quantile mapping techniques in the multivariate case. We also propose a definition of non-stationary bias correction as a transfer of the model to the observational world, and we extend our method in this context. Those methodologies are first tested on an idealized chaotic system with three variables. In those controlled experiments, the correlations between variables appear almost perfectly corrected by our method, as opposed to a univariate correction. Our methodology is also tested on daily precipitation and temperatures over 12 locations in southern France. The correction of the inter-variable and inter-site structures of temperatures and precipitation appears in agreement with the multi-dimensional evolution of the model, hence satisfying our suggested definition of non-stationarity.