Currently I am developing an online algorithm for harmonic tidal prediction. The measurements collected are $(y_i , t_i)$, with $y_i$ the height of the water (in meter) and $t_i$ the time of measurement (in hours). As the tide is (assumed to be) a result of $d$ harmonic effects, it can be modelled as
$$ y \approx \sum_{j=1} ^{d} w^j \cos(H_j t ) + \sum_{j = 1} ^{d} w^{d+j} \sin(H_j t) $$ Where $(H_1 \cdots H_d)$ are constants corresponding with the harmonic effects, and $(w^1 \cdots w^{2d} )$ the weights to be determined. If we define $x_i$ as $$x_i = \begin{pmatrix} \cos(H_1 t_i)\\ \vdots \\ \cos(H_d t_i) \\ \sin(H_1 t_i) \\ \vdots \\ \sin(H_d t_i) \end{pmatrix} $$
On this transformed data $(y_i, x_i)$ I employ a recursive least squares algorithm, characterised by the following formulae: $$ \Gamma_i = \Gamma_{i-1} - \dfrac{\Gamma_{i-1}x_i x_i ^T \Gamma_{i-1}}{1+ x_i ^T\Gamma_{i-1} x_i} \\ w_i = w_{i-1} - \Gamma_{i}x_i (x_i ^T w_{i-1} -y_i) $$ Where $w_0 = 0 \in \mathbb{R}^{2d}$ and $\Gamma_0 = 0 \in \mathbb{R}^{2d \times 2d}$.
I have the algorithm working for $d \leq 5$, yet at a try for $d = 37$, a flaw has appeared.
The predictions get off to a wrong start, predicting high when low water and vice versa. As the algorithm's correction of $w_i$ is proportionate to the prediction error (due to the last factor in the definition of $w_i$), it starts to overcompensate, resulting in an even bigger error in the next step. Then this big error induces another big overcompensation in the next step, etc. etc...
A little illustration is shown below, illustrating the abnormal growth in predictions between steps 1, 2 and 3.
My question is if someone has any thought on what the cause might be and ,if possible, what could be a solution to this.
PS: I realize this is a very specific question, yet I hope I have given sufficient information on the problem without stretching this post any further. If not, please comment.