De forma similar a la respuesta de Juan Gil pero un poco diferente yo diría lo siguiente basándome en este :
El proceso de la OU $$dX_t = \kappa(\theta-X_t)dt + \sigma dW_t$$ se puede discretizar (discretización de Euler-Maryuama) en tiempos $n \Delta t,n=1,\ldots,\infty $ que da con $t = k \Delta t$ $$ X_{k+1} - X_k = \kappa \theta \Delta t -\kappa X_k \Delta t + \sigma (W_{k+1} - W_k), $$ reordenación y ajuste $\sigma (W_{k+1} - W_k) = \sigma \sqrt{\Delta t} \epsilon_k $ nos encontramos con que: $$ X_{k+1} = \kappa \theta \Delta t - (\kappa \Delta t - 1) X_k + \sigma \sqrt{\Delta t} \epsilon_k. $$ Así que puede modelar un proceso AR(1) y luego identificar los parámetros utilizando la ecuación anterior.
Pensando en ello de nuevo uno puede probablemente dejar $X_{k+1} - X_k$ en la lhs y entonces uno simplemente hace una regresión pero no sé exactamente sobre los términos de error en este caso.
He encontrado esto con Código R, allí se utiliza un enfoque MLE. Encontrará varias soluciones en este Stack Overflow pregunta .