Since you are talking about phase difference between two time series, say $\{x_k\}$and$\{y_k\}$ the underlying signal must be something that can be viewed as a "sinusoid" in some approximate sense. If you know nothing else then with a Fourier transform you can estimate the underlying sinusoidal frequency, say you get $\omega$. Now apply to the time series what RF engineers call an IQ (coherent) detection: calculate two sums of products:$ M(x) = \sum_k x_k e^{-\mathfrak j \omega t_k}$ and $ M(y) = \sum_k y_k e^{-\mathfrak j \omega t_k}$, and let $z=M(x) \bar M(y)$ then the phase of $z$ is the estimate of the phase angle between the two time series.
This can be seen as follows. Let $r_k = \rm{cos}(\omega t_k + \phi)+\nu_k$ where $\nu_k$ is noise.
Then $$ M(r) = \sum_k \rm{cos}(\omega t_k + \phi_r) e^{-\mathfrak j \omega t_k} +\rm{noise}\\
=\frac{1}{2}\sum_k [e^{\mathfrak j (\omega t_k + \phi_r)}+e^{-\mathfrak j (\omega t_k + \phi)}] e^{-\mathfrak j \omega t_k}+\rm{noise}\\
=\frac{1}{2}Ne^{\mathfrak j \phi_r} + \sum_k e^{-\mathfrak j (2\omega t_k + \phi_r)}+\rm{noise}\\$$
The sum is nearly zero if many samples are taken over a period, and it is exactly zero if the period is exactly sampled. The rms of the independent noise samples grows like $\sqrt{N}$. In this context $\phi$ is the phase of the signal $r$ relative to the receiver's oscillator. So if you calculate the phases for both signals relative to the local clock (oscillator, sample) then the differecne of those will get you the relative phase between them:$M(x)\bar M(y) \propto e^{\mathfrak j (\phi_x - \phi_y)}$
