  We propose Multivariate Quantile Function Forecaster (MQF$^2$), a global
probabilistic forecasting method constructed using a multivariate quantile
function and investigate its application to multi-horizon forecasting. Prior
approaches are either autoregressive, implicitly capturing the dependency
structure across time but exhibiting error accumulation with increasing
forecast horizons, or multi-horizon sequence-to-sequence models, which do not
exhibit error accumulation, but also do typically not model the dependency
structure across time steps. MQF$^2$ combines the benefits of both approaches,
by directly making predictions in the form of a multivariate quantile function,
defined as the gradient of a convex function which we parametrize using
input-convex neural networks. By design, the quantile function is monotone with
respect to the input quantile levels and hence avoids quantile crossing. We
provide two options to train MQF$^2$: with energy score or with maximum
likelihood. Experimental results on real-world and synthetic datasets show that
our model has comparable performance with state-of-the-art methods in terms of
single time step metrics while capturing the time dependency structure.
