  We study the pulsar timing properties and the data analysis methods during
glitch recoveries. In some cases one first fits the time-of-arrivals (TOAs) to
obtain the "time-averaged" frequency $\nu$ and its first derivative $\dot\nu$,
and then fits models to them. However, our simulations show that $\nu$ and
$\dot\nu$ obtained this way are systematically biased, unless the time
intervals between the nearby data points of TOAs are smaller than about $10^4$
s, which is much shorter than typical observation intervals. Alternatively,
glitch parameters can be obtained by fitting the phases directly with
relatively smaller biases; but the initial recovery timescale is usually chosen
by eyes, which may introduce a strong bias. We also construct a
phenomenological model by assuming a pulsar's spin-down law of $\dot\nu\nu^{-3}
=-H_0 G(t)$ with $G(t)=1+\kappa e^{-t/\tau}$ for a glitch recovery, where $H_0$
is a constant and $\kappa$ and $\tau$ are the glitch parameters to be found.
This model can reproduce the observed data of slow glitches from B1822--09 and
a giant classical glitch of B2334+61, with $\kappa<0$ or $\kappa>0$,
respectively. We then use this model to simulate TOA data and test several
fitting procedures for a glitch recovery. The best procedure is: 1) use a very
high order polynomial (e.g. to 50th order) to precisely describe the phase; 2)
then obtain $\nu(t)$ and $\dot\nu(t)$ from the polynomial; and 3) the glitch
parameters are obtained from $\nu(t)$ or $\dot\nu(t)$. Finally, the uncertainty
in the starting time $t_0$ of a classical glitch causes uncertainties to some
glitch parameters, but less so to a slow glitch and $t_0$ of which can be
determined from data.
