We study the pulsar timing properties and the data analysis methods of glitch recoveries. In some cases one first fits the times of arrival (TOAs) to obtain the "time-averaged" frequency. and its first derivative.., and then fits models to them. However, our simulations show that. and. obtained in 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 time scale is usually chosen by eye, which may introduce a strong bias. We also construct a phenomenological model by assuming a pulsar spin-down law of (nu) over dot nu(-3) = -H(0)G(t) with G(t) = 1 + kappa e(-t/iota) for a glitch recovery, where H-0 is a constant and. and t are the glitch parameters to be found. This model can reproduce the observed data of slow glitches from B1822-09 and the 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) to use a very high order polynomial (e. g., to 50th order) to precisely describe the phase, (2) obtain.(t) and..(t) from the polynomial, then (3) obtain the glitch parameters from nu(t) or (nu) over dot (t). Finally, the uncertainty in the starting time t(0) of a classical glitch causes uncertainties in some glitch parameters, but less so for a slow glitch, t(0), of which can be determined from data.