A novel approximate technique based on Wiener path integrals (WPIs) is developed for determining, in a computationally efficient manner, the non-stationary joint response probability density function (PDF) of nonlinear multi-degree-of-freedom dynamical systems. Specifically, appropriate multi-dimensional (time-dependent) global bases are constructed for approximating the non-stationary joint response PDF. In this regard, two distinct approaches are pursued. The first employs expansions based on Kronecker products of bases (e.g., wavelets), while the second utilizes representations based on positive definite functions. Next, the localization capabilities of the WPI technique are exploited for determining PDF points in the joint space-time domain to be used for evaluating the expansion coefficients at a relatively low computational cost. In contrast to earlier implementations of the WPI technique, the herein developed generalization and enhancement circumvents computationally daunting brute-force discretizations of the time domain in cases where the objective is to determine the complete time-dependent non-stationary response PDF. Several numerical examples pertaining to diverse structural systems are considered, including both single- and multi-degree-of-freedom nonlinear dynamical systems subject to non-stationary excitations, as well as a bending beam with a non-Gaussian and non-homogeneous Young's modulus. Comparisons with pertinent Monte Carlo simulation data demonstrate the accuracy of the technique. (C) 2019 Elsevier Ltd. All rights reserved.