  This paper presents rigorous tests of pulsar timing array methods and
software, examining their consistency across a wide range of injected
parameters and signal strength. We discuss updates to the 15-year isotropic
gravitational-wave background analyses and their corresponding code
representations. Descriptions of the internal structure of the flagship
algorithms \texttt{Enterprise} and \texttt{PTMCMCSampler} are given to
facilitate understanding of the PTA likelihood structure, how models are built,
and what methods are currently used in sampling the high-dimensional PTA
parameter space. We introduce a novel version of the PTA likelihood that uses a
two-step marginalization procedure that performs much faster when the white
noise parameters remain fixed. We perform stringent tests of consistency and
correctness of the Bayesian and frequentist analysis software. For the Bayesian
analysis, we test prior recovery, injection recovery, and Bayes factors. For
the frequentist analysis, we test that the cross-correlation-based optimal
statistic, when modified to account for a non-negligible gravitational-wave
background, accurately recovers the amplitude of the background. We also
summarize recent advances and tests performed on the optimal statistic in the
literature from both GWB detection and parameter estimation perspectives. The
tests presented here validate current and future analyses of PTA data.
