We propose a methodology to characterize a multivariate non-stationary vector random process that can be used for simulating random realizations that keep the probabilistic behavior of the original time series. The marginal probability distribution of each component process is assumed to be a piecewise function deﬁned by several weighted parametric probability models. The weights are obtained analytically by ensuring that the probability density function is well deﬁned and that it is continuous at the common endpoints. The probability model is assumed to vary periodically in time over a predeﬁned time period by deﬁning the model parameters and the common endpoints as truncated generalized Fourier series. The coeﬃcients of the expansions are obtained with the maximum likelihood method. Three diﬀerent types of sets of orthogonal functions are tested. The method is applied to three time series with diﬀerent particularities. Firstly, it is shown its good behavior to capture the highly variable freshwater discharges at a dam located in a semiarid zone in Andalucía (Spain) which is inﬂuenced not only by the climate variability but also by management decisions. Secondly, for the Wolf sunspot number time series, the Schwabe cycle and time variations close to the 7.5 and 17 years are analyzed along a 22-year cycle. Finally, the method is applied to a bivariate (velocity and direction) wind time series observed at a location of the Atlantic Ocean. For this case, the analysis, that was combined with a vectorial autoregresive model, focus on the assessment of the goodness of the methodology to replicate the statistical features of the original series. In particular, it is found that it reproduces the marginal and joint distributions, the wind rose, and the duration of sojourns above given thresholds.