  Although it is widely understood that pulsar timing observations generally
contain time-correlated stochastic signals (TCSSs; red timing noise is of this
type), most data analysis techniques that have been developed make an
assumption that the stochastic uncertainties in the data are uncorrelated, i.e.
"white". Recent work has pointed out that this can introduce severe bias in
determination of timing-model parameters, and that better analysis methods
should be used. This paper presents a detailed investigation of timing-model
fitting in the presence of TCSSs, and gives closed expressions for the post-fit
signals in the data. This results in a Bayesian technique to obtain
timing-model parameter estimates in the presence of TCSSs, as well as
computationally more efficient expressions of their marginalised posterior
distribution. A new method to analyse hundreds of mock dataset realisations
simultaneously without significant computational overhead is presented, as well
as a statistically rigorous method to check the internal consistency of the
results. As a by-product of the analysis, closed expressions of the rms
introduced by a stochastic background of gravitational-waves in
timing-residuals are obtained. Using $T$ as the length of the dataset, and
$h_c(1\rm{yr}^{-1})$ as the characteristic strain, this is: $\sigma_{\rm GWB}^2
= h_{c}(1\rm{yr}^{-1})^2 (9\sqrt[3]{2\pi^4}\Gamma(-10/3) / 8008) \rm{yr}^{-4/3}
T^{10/3}$.
