Similar to Juan Gil's answer but a bit differently I would say the following based on this:
The OU process $$dX_t = \kappa(\theta-X_t)dt + \sigma dW_t$$ can be (Euler-Maryuama discretization) discretized at times $n \Delta t,n=1,\ldots,\infty $ which gives with $t = k \Delta t$
$$
X_{k+1} - X_k = \kappa \theta \Delta t -\kappa X_k \Delta t + \sigma (W_{k+1} - W_k),
$$
rearranging and setting $\sigma (W_{k+1} - W_k) = \sigma \sqrt{\Delta t} \epsilon_k $ we get:
$$
X_{k+1} = \kappa \theta \Delta t - (\kappa \Delta t - 1) X_k + \sigma \sqrt{\Delta t} \epsilon_k.
$$
So you can model an AR(1) process and then identify the parameters using the equation above.
Thinking about it again one can probably leave
$X_{k+1} - X_k$ on the lhs and then one simply does a regression but I don't know exactly about the error terms in this case.
I have found this with R code, there an MLE approach is used. You find various solutions in this Stack Overflow question.