There are two ways this is commonly handled, both still based on O-U:
- Deterministic long-term fluctuations.
- Separate timescale processes
Let's say your real time series of interest is $Y(t)$. In the first case, often used by econometricians to handle seasonality, one sets up a deterministic process $F(t)$, and then run Ornstein-Uhlenbeck as a variation on that process. That is to say, we define the stochastic variable $X(t)=Y(t)-F(t)$ and then specify
$$
dX_t = -\alpha_t X_t + \sigma_t dW
$$
The second case is more physical, and one uses a two-dimensional OU process
$$
dY^{(1)}_t = -\alpha_t^{(1)} (Y^{(1)}_t - \bar{Y})+ \sigma^{(1)}_t dW_1
$$
$$
dY^{(2)}_t = -\alpha_t^{(2)} Y^{(2)}_t + \sigma^{(2)}_t dW_2
$$
with $\alpha_1 \ll \alpha_2$, after which we define $Y=Y^{(1)}+Y^{(2)}$.
The math is all still quite tractable, particularly so with constant $\alpha,\sigma$.
It's worth noting that, with respect to calibration, you can treat the problem as separable to the extent that $\alpha_1 \ll \alpha_2$ is really true.
That is to say, you get a far better fit using a Fourier transform. Use a low-pass filter to obtain data for finding $\sigma^{(1)}$ and $\alpha^{(1)}$. Then, you use a high-pass filter to get $\sigma^{(2)}$ and $\alpha^{(2)}$.
For example assume $\alpha^{(1)}, \sigma^{(1)}$ are roughly daily and $\alpha^{(2)}, \sigma^{(2)}$ are roughly annual. We define
$$
f(\omega)={\cal{F}}(Y)
$$
and let
$$
f^{(1)}(\omega) = f(\omega) \times \mathbb{1} \lbrace \omega<(30d)^{-1} \rbrace
$$
and
$$
f^{(2)}(\omega) = f(\omega) \times \mathbb{1} \lbrace \omega\geq(30d)^{-1} \rbrace,
$$
taking due care of the degeneracies in fourier transforms for real data sets. Now set
$$
X^{(1)} = \text{Re} \left[ \cal{F}^{-1}\left( f^{(1)} \right) \right]
$$
and
$$
X^{(2)} = \text{Re} \left[ \cal{F}^{-1}\left( f^{(2)} \right) \right].
$$
Running a standard maximum-likelihood estimate on $X^{(2)}$ will find $\sigma^{(2)}$ and $\alpha^{(2)}$. Subsampling $X^{(2)}$ on, say, the same 30 day interval as the cutoff will give a time series appropriate for estimating $\sigma^{(1)}$ and $\alpha^{(1)}$.
Running these two separate two-dimensional calibrations (say as 2-dimensional maximum-likelihood estimates) is far more stable than a single four-dimensional calibration. You still have to be careful. For example if $\sigma^{(2)}$ is too small you have no hope of finding it in all the noise of the lower-frequency process.