We present a 5-step method to estimate piecewise linear trends in SBUV MOD v6 monthly zonal mean (MZM) profile ozone timeseries containing gaps (e.g. associated with volcanic perturbations). The method is non-parametric and does not use proxies. In Step 1, an optimised local regression (LOESS) to the porous data is obtained using bootstrap in a way that guarantees trend stationarity via the KPSS hypothesis test. In Step 2, the discrete Fourier transform is computed for the longest gap-free segment of the LOESS-centered data to identify the highest relative variance frequencies in the data. In Step 3, these Fourier peaks are used together with the LOESS trend to construct a sum of sines fit to the full timeseries and interpolate any missing data. In Step 4, Monte-Carlo Singular Spectrum Analysis (MC-SSA) is applied to the reconstructed timeseries to calculate a (5,95 percent) confidence interval around the median long-term nonlinear trend obtained from a bootstrap sample of 2000 runs. In Step 5, a turnaround point is set at January 1997 and used to regress a 2-segment PWLT fit pre- and post-1997 to the MC-SSA trend using least squares optimisation. We present data-driven profile ozone trend estimates (percent / decade) from 25.45 hPa to the top of atmosphere in the equatorial bands (20S-20N and 50S-50N) and in mid-latitude bands (50S-35S and 35N-50N) for comparison with proxy-driven approaches. The maximum negative trend detected is -7.41 percent / decade (pre-1997) in the Northern mid-latitude band in the 2.55 - 1.61 hPa pressure layer which switches to a positive trend of 1.74 percent / decade post-1997. In addition to providing predicted values for missing months in the dataset, Fourier peak frequencies extracted can also be used to test for the presence of known cyclical proxies. Some results from wavelet analysis are presented to illustrate.