While we are usually focused on forecasting future values of time series, itis often valuable to additionally predict their entire probabilitydistributions, e.g. to evaluate risk, Monte Carlo simulations. On example oftime series of $\approx$ 30000 Dow Jones Industrial Averages, there will bepresented application of hierarchical correlation reconstruction for thispurpose: mean-square estimating polynomial as joint density for (current value,context), where context is for example a few previous values. Then substitutingthe currently observed context and normalizing density to 1, we get predictedprobability distribution for the current value. In contrast to standard machinelearning approaches like neural networks, optimal polynomial coefficients herecan be inexpensively directly calculated, have controllable accuracy, areunique and independent, each has a specific cumulant-like interpretation, andsuch approximation using can approach complete description of any real jointdistribution - providing a perfect tool to quantitatively describe and exploitstatistical dependencies in time series. There is also discussed applicationfor non-stationary time series: adapting coefficients to local statisticalbehavior.