Supermassive black hole binaries (SMBHBs) are expected to emit continuous gravitational waves in the pulsar timing array (PTA) frequency band ($10^{-9}$–$10^{-7}$ Hz). The development of data analysis techniques aimed at efficient detection and characterization of these signals is critical to the gravitational wave detection effort. In this paper we leverage methods developed for LIGO continuous wave gravitational searches, and explore the use of the $\mathcal{F}$-statistic for such searches in pulsar timing data. Babak & Sesana 2012 have already used this approach in the context of PTAs to show that one can resolve multiple SMBHB sources in the sky. Our work improves on several aspects of prior continuous wave search methods developed for PTA data analysis. The algorithm is implemented fully in the time domain, which naturally deals with the irregular sampling typical of PTA data and avoids spectral leakage problems associated with frequency domain methods. We take into account the fitting of the timing model, and have generalized our approach to deal with both correlated and uncorrelated colored noise sources. We also develop an incoherent detection statistic that maximizes over all pulsar dependent contributions to the likelihood. To test the effectiveness and sensitivity of our detection statistics, we perform a number of monte-carlo simulations. We produce sensitivity curves for PTAs of various configurations, and outline an implementation of a fully functional data analysis pipeline. Finally, we present a derivation of the likelihood maximized over the gravitational wave phases at the pulsar locations, which results in a vast reduction of the search parameter space.