  Pulsar timing arrays (PTAs) detect gravitational waves (GWs) via the
correlations they induce in the arrival times of pulses from different pulsars.
We assume that the GWs are described by a Gaussian ensemble. The mean
correlation $h^2 \mu_{\rm u}(\gamma)$ as a function of the angle $\gamma$
between the directions to two pulsars was predicted by Hellings and Downs (HD)
in 1983. The variance $\sigma^2_{\rm tot}(\gamma)$ in this correlation was
recently calculated by Allen[11] for a single noise-free pulsar pair at angle
$\gamma$, which shows that after averaging over many pairs, the variance
reduces to an intrinsic cosmic variance $\sigma^2_{\rm cos}(\gamma)$. Here, we
extend this to an $arbitrary$ set of pulsars at specific sky locations, with
pulsar pairs binned by $\gamma$. We derive the linear combination of
pulsar-pair correlations which is the optimal estimator of the HD correlation
for each bin, illustrating our methods with plots of the expected range of
variation away from the HD curve, for the sets of pulsars monitored by three
active PTA collaborations. We compute the variance of and the covariance
between these binned estimates, and show that these reduce to the cosmic
variance and covariance $s(\gamma,\gamma')$ respectively, in the many-pulsar
limit. The likely fluctuations away from the HD curve $\mu_{\rm u}(\gamma)$ are
strongly correlated/anticorrelated in the three angular regions where $\mu_{\rm
u}(\gamma)$ is successively positive, negative, and positive. We also construct
the optimal estimator of the squared strain $h^2$. When there are very many
pulsar pairs, this determines $h^2$ with arbitrary precision because PTAs probe
an infinite set of GW modes. To assess observed deviations away from the HD
curve, we characterize several $\chi^2$ goodness-of-fit statistics. We also
show how pulsar noise and measurement noise can be included.
