What is the explanation for the inverse FFT of a cross-correlated signal being shifted by some amount in the time domain?
The issue is that the FFT version always assumes that the indices start numbering from $n=0$.  The result of the FFT version will be periodic, with period the same as the FFT length.

To show the same result as the correlation version, you need to use numpy.fft.fftshift.

Here is my attempt to generate your original signals:



Then the two versions of correlation:



which seems to do what you're seeing.

Whnn I apply fftshift to the FFT version, I get:





Code Only Below

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
t = np.linspace(-0.2, 0.5, 500, endpoint=False)
t2 = np.linspace(-0.2, 0.5, 250, endpoint=False)
t3 = np.linspace(-0.4, 1.0, 999)
square = (signal.square(2 * np.pi * t - np.pi/2 , duty=0.2) + 1)/3
triangle = np.zeros(500)
triangle[250:500] = (signal.sawtooth(2 * np.pi * t2 - np.pi, 0 ) + 1)/2

txcorr = signal.correlate(square, triangle, mode = 'full', method='direct')
txcorr2 = np.fft.ifft(np.fft.fft(square, 999)*np.conj(np.fft.fft(triangle, 999)))

plt.figure(0)
plt.plot(t, square)
plt.plot(t, triangle)
plt.figure(1)
plt.plot(t3, txcorr)
plt.plot(t3, txcorr2)
plt.figure(2)
plt.plot(t3, np.fft.fftshift(txcorr2), '.')
plt.plot(t3, txcorr)