  We model the evolution of the spin frequency's second derivative $\ddot\nu$
and the braking index $n$ of radio pulsars with simulations within the
phenomenological model of their surface magnetic field evolution, which
contains a long-term power-law decay modulated by short-term oscillations. For
the pulsar PSR B0329+54, a model with three oscillation components can
reproduce its $\ddot\nu$ variation. We show that the "averaged" $n$ is
different from the instantaneous $n$, and its oscillation magnitude decreases
abruptly as the time span increases, due to the "averaging" effect. The
simulated timing residuals agree with the main features of the reported data.
Our model predicts that the averaged $\ddot\nu$ of PSR B0329+54 will start to
decrease rapidly with newer data beyond those used in Hobbs et al.. We further
perform Monte Carlo simulations for the distribution of the reported data in
$|\ddot\nu|$ and $|n|$ versus characteristic age $\tau_{\rm c}$ diagrams. It is
found that the magnetic field oscillation model with decay index $\alpha=0$ can
reproduce the distributions quite well. Compared with magnetic field decay due
to the ambipolar diffusion ($\alpha=0.5$) and the Hall cascade ($\alpha=1.0$),
the model with no long term decay ($\alpha=0$) is clearly preferred for old
pulsars by the p-values of the two-dimensional Kolmogorov-Smirnov test.
