While we are usually focused on forecasting future values of time series, it is often valuable to additionally predict their entire probability distributions, e.g. to evaluate risk, Monte Carlo simulations. On example of time series of $\approx$ 30000 Dow Jones Industrial Averages, there will be presented application of hierarchical correlation reconstruction for this purpose: mean-square estimating polynomial as joint density for (current value, context), where context is for example a few previous values. Then substituting the currently observed context and normalizing density to 1, we get predicted probability distribution for the current value. In contrast to standard machine learning approaches like neural networks, optimal polynomial coefficients here can be inexpensively directly calculated, have controllable accuracy, are unique and independent, each has a specific cumulant-like interpretation, and such approximation using can approach complete description of any real joint distribution - providing a perfect tool to quantitatively describe and exploit statistical dependencies in time series. There is also discussed application for non-stationary time series like calculating linear time trend, or adapting coefficients to local statistical behavior.