Adding to @Dilip Sarwate's dot product solution, which gives you the magnitude of the phase difference (note that arccos
returns a value from [0, $\pi$]), if you want to know which signal leads the other, you'll need an expression that includes $sin$, such as the cross product.
$$\left\|\mathbf{a} \times \mathbf{b}\right\| = \left\| \mathbf{a} \right\| \left\| \mathbf{b} \right\| \sin (\theta)$$
$$\sin \theta = \frac{ \left\|\mathbf{a} \times \mathbf{b}\right\| }{\left\| \mathbf{a} \right\| \cdot \left\| \mathbf{b} \right\|} = \frac{\displaystyle \sum_{n=0}^{N-1} t_x[n]*y[n] - t_y[n]*x[n]}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}(x[n])^2\displaystyle \sum_{n=0}^{N-1}(y[n])^2}}$$
with $\mathbf{a} = \left< t, x,0\right>$ and $\mathbf{a} = \left< t, y,0\right>$
Which when using the same time base for the 2 signals, simplifies to:
$$\sin \theta = \frac{\displaystyle \sum_{n=0}^{N-1} t[n]\left(y[n] - x[n]\right)}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}(x[n])^2\displaystyle \sum_{n=0}^{N-1}(y[n])^2}}$$
In Python, without knowing the phase shift magnitude or direction in advance:
t = np.deg2rad(np.arange(0, 360*4))
y1 = np.sin(t)
y2 = np.sin(t+10*np.pi/180)
# from dot product
opp = np.sum(y2*y1)
hyp = np.sqrt( np.sum(y2**2) * np.sum(y1**2) )
# from cross product
adj = np.cross([np.cos(t), y1], [np.cos(t), y2], axis=0)
phase_angle = np.atan2( adj, opp )
Which is consistent with $\sin(\theta) = opp/hyp$, $\cos(\theta) = adj/hyp$, and $\tan(\theta) = opp/adj$,