Thank you for your interest. The identification algorithm uses this solution:
Basically consists in integrating the inputs and outputs of the system (according to the desired model order) and then solving a couple of (linear) least squares problems.
The delay in nonlinear, so given that the (linear) least squares problems are “relatively cheap”, the delay is calculated by brute force initially using a grid of values. Then refinement is made using a simple newton method.
Hope this is helpful.