Pulsar Timing Array (PTA) experiments aim to detect gravitational wave (GW) signals at nanohertz frequencies through pulsar timing, by characterising the timing delays they are expected to induce. Pulsars are incredibly stable, compact, fast-rotating objects with a collimated jet of radio emission, observed as a pulsating signal (lighthouse effect). GWs are expected to induce spatially correlated delays in an array of pulsar observations. This cross-pulsar correlation is modelled by the Hellings and Downs (HD) curve. The most likely signal in the PTA band is a stochastic GW background (GWB), due to the incoherent superposition of monochromatic GW emissions from supermassive black hole binaries (SMBHBs). This thesis consists of several works describing the impact of realistic and astrophysically motivated GW source models on GWB signal recovery, using both analytical predictions and simulation studies, and the development of new methods to both improve the efficiency of current PTA data analysis pipelines, and finally extend PTA experiments to gamma-ray pulsars. Firstly, the variance in the HD recovery is quantified analytically for different cases of finite ensembles of discrete binary GW sources, with an aim to investigate the impact of various effects: (i) the discreteness of the sources in the ensemble; (ii) the arbitrary inclination of the GW source orbit with respect to the line of sight; (iii) the presence of correlations in their sky locations. While the obtained mean correlation is always proportional to the HD curve, the variance increases with the degrees of freedom added to the GW source models. However, in the limit of infinite sources, the variance converges to the variance obtained for a simple Gaussian ensemble of GW sources. The following chapter then presents a novel regularized likelihood formulation in Fourier space. This alternative method has many advantages. It divides the GWB search into a two-step analysis, the first step of which can be carried out in parallel for all pulsars individually. Furthermore, it can also be applied to gamma-ray photon data and paves the way for searches of cross-pulsar correlated signals from combined radio and gamma-ray PTA data. This method retains all the information regarding the signals of interest, eliminating the complications of underlying noise models and implementation differences. The limitations, strengths, potential, and possible intrinsic biases of the current PTA data analysis methods were then tested in detail using simulation studies. Particular focus was given to (a) studying the impact of the assumption that the GWB spectra is well modelled by a simple power law (which can prevent us from correctly characterising the SMBHB population that produced the signal); (b) developing and testing a new method to distinguish between circular and eccentric SMBHB populations, by quantifying the correlation between skymaps of GW power at different frequency. Finally, a novel, optimised algorithm for timing gamma-ray pulsars following a Gibbs sampling scheme is presented and tested for statistical robustness. By applying this method and the regularized likelihood formulation to the first gamma-ray PTA dataset, it was possible to carry out the first search for an HD correlated signal directly on pulsars' gamma-ray photon data and derive an updated and more robust upper limit on the amplitude of a GWB signal for the gamma-ray PTA.