An Ornstein-Uhlenbeck (OU) process can be considered as a continuous time interpolation of the discrete time AR$(1)$ process. Departing from this fact, we analyse in this work the effect of iterating OU treated as a linear operator that maps a Wiener process onto Ornstein-Uhlenbeck process, so as to build a family of higher order Ornstein-Uhlenbeck processes, OU$(p)$, in a similar spirit as the higher order autoregressive processes AR$(p)$. We show that for $p \ge 2$ we obtain in general a process with covariances different than those of an AR$(p)$, and that for various continuous time processes, sampled from real data at equally spaced time instants, the OU$(p)$ model outperforms the appropriate AR$(p)$ model. Technically our composition of the OU operator is easy to manipulate and its parameters can be computed efficiently because, as we show, the iteration of OU operators leads to a process that can be expressed as a linear combination of basic OU processes. Using this expression we obtain a closed formula for the covariance of the iterated OU process, and consequently estimate the parameters of an OU$(p)$ process by maximum likelihood or, as an alternative, by matching correlations, the latter being a procedure resembling the method of moments.