I am not sure to understand why you would like to use BIC, as it's usually a tool for model selection rather than a meaningful statistics for timeseries.

Another approach that works quite well could be to smooth your noisy signal to remove the noise (via a mooving average for instance), and remove the trend in your signal. Then use a Fourier transform and/or correlation to detect the periodicity in the signal (which should be the period at which the mean changes). From this it should be easy to approximate the means.

Here is a small example that I tested which worked quite well as a first approximation:

from scipy.fft import fft, fftfreq
n = 50
ma = np.convolve(y,np.ones(n), mode='valid')/n # denoised signal
rm_trend = y-((ma[-1]-ma[0])/len(ma)*np.arange(len(y))+ma[0]) # remove trend

corr = np.correlate(rm_trend,rm_trend,mode='full')
corr = corr[corr.shape[0]//2:]

y_fft = fft(y,norm='forward')[1:len(y)//2] # remove the mean

corr = np.correlate(rm_trend,rm_trend,mode='full') # autocorrelation
corr = corr[corr.shape[0]//2:]
freq = fftfreq(len(corr)) # frequencies
corr_fft = fft(corr,norm='forward')[1:len(corr)//2] # FFT without mean
k = 1/freq[np.argmax(corr_fft )+1]
print(k)

Please tell me if this does not answer your question.