I added the following comment to the code so folks can better understand what is going on here. You can see the original Levinson-Durbin Recursion algorithm in the commented out note.
// Original Levinson-Durbin algorithm used to implement Levinson recursion
// where a - coefficients of the model, p - order of the model.
// Here we need to find the autoregressive coefficients by solving directly
// our set of equations with n=2*p by the Levinson-Durbin method. Such method
// of prediction is called Prony Method; however, its disadvantage is the
// instability during the prediction of the future values of the series. That's
// why this method has not been included and instead we use a modified
// Levinson Recursion to calculate the prediction coefficients.
// I've included the origina method so one can compare the differences. You'll
// notice that both methods are very similar but the modified version gives the
// desired results. The difference is that the modified version calculates the
// coefficients a[] by decreasing the mean-root-square error on the training
// last n-p bars