We model the evolution of the magnetic fields of neutron stars as consisting of a long term power-law decay modulated by short term small amplitude oscillations. Our model predictions on the timing noise of neutron stars agree well with the observed statistical properties and correlations of normal radio pulsars. For individual pulsars our model can effectively reduce their timing residuals, thus offering the potential of more sensitive detections of gravitational waves with pulsar timing arrays. Finally our model can also re-produce their observed correlation and oscillations of second derivative of frenquency, as well as the “slow glitch” phenomenon.