  A new Bayesian method for the analysis of folded pulsar timing data is
presented that allows for the simultaneous evaluation of evolution in the pulse
profile in either frequency or time, along with the timing model and additional
stochastic processes such as red spin noise, or dispersion measure variations.
We model the pulse profiles using `shapelets' - a complete ortho-normal set of
basis functions that allow us to recreate any physical profile shape. Any
evolution in the profiles can then be described as either an arbitrary number
of independent profiles, or using some functional form. We perform simulations
to compare this approach with established methods for pulsar timing analysis,
and to demonstrate model selection between different evolutionary scenarios
using the Bayesian evidence. %s The simplicity of our method allows for many
possible extensions, such as including models for correlated noise in the pulse
profile, or broadening of the pulse profiles due to scattering. As such, while
it is a marked departure from standard pulsar timing analysis methods, it has
clear applications for both new and current datasets, such as those from the
European Pulsar Timing Array (EPTA) and International Pulsar Timing Array
(IPTA).
