I also think that the diff equations are better suited if you break the time into two component of {year, month}.

In the .mml equations I gave from the SVR's Gaussian Wavelet you see {x1, x2} accordingly.

You do not need to make any mod 12 or discrete year, assume two continuous vars.

My hunch is that the equations produce much better models/solutions if the time is broken into 1 period value and another linear.

You could also add the leap year as a param by itself 1 or 0, and maybe other params like # of sunspots and so on: {x1, x2, x3 ...} and I could compute you a closed form equation from SVR regressor just the same.