Geodetic time series are usually modeled with a deterministic approach that includes trend, annual, and semiannual periodic components having constant amplitude and phase-lag. Although simple, this approach neglects the time-variability or stochasticity of trend and seasonal components, and can potentially lead to inadequate interpretations, such as an overestimation of global navigation satellite system (GNSS) station velocity uncertainties, up to masking important geophysical phenomena. In this contribution, we generalize previous methods for determining trends and seasonal components and address the challenge of their time-variability by proposing a novel linear additive model, according to which (i) the trend is allowed to evolve over time, (ii) the seasonality is represented by a fractional sinusoidal waveform process (fSWp), accounting for possible non-stationary cyclical long-memory, and (iii) an additional serially correlated noise captures the short term variability. The model has a state space representation, opening the way for the evaluation of the likelihood and signal extraction with the support of the Kalman filter (KF) and the associated smoothing algorithm. Suitable enhancements of the basic methodology enable handling data gaps, outliers, and offsets. We demonstrate the advantage of our method with respect to the benchmark deterministic approach using both observed and simulated time series and provide a fair comparison with the Hector software. To that end, various geodetic time series are considered which illustrate the ability to capture the time-varying stochastic seasonal signals with the fSWp.