Why is the bias instability read from my Allan Deviation plot different than my simulation input? [closed]

Question

I have simulated a gyroscope at rest with angular random walk and bias instability, following the mathworks IMU simulation model found here. In conjunction with that simulation model, I am plotting the Allan Deviation and trying to derive the input parameters for random walk and instability from the plot following another mathworks resource.

For the case of the ARW calculated from Allan deviation, I'm getting back my input parameters which is great. However, when it comes to the instability term, my input and the Allan deviation is off by an order of magnitude, and I can't figure out where I went wrong.

D2R = np.pi / 180
fs = 1200               # frequency - Hz
T = 1 / fs              # period - sec
runtime = 3600          # sec
num_samples = fs * runtime

x = np.arange(0, runtime, 1 / fs)  # seconds
omega = np.zeros_like(x)

#####################################
#### SIMULATE ANGULAR RANDOM WALK ###
#####################################
arw = 0.05                        # deg/rt-hour
arw = arw * D2R / 60              # rad/rt-second
arw = arw * np.sqrt(fs)           # rad/sec
arw = arw * np.random.randn(num_samples)


##################################
#### SIMULATE BIAS INSTABILITY ###
##################################
gyro_bias_stability = 4            # deg/hour
gyro_bias_stability *= D2R / 3600  # rad/sec

ABS = np.zeros_like(omega)
for i in range(1, len(ABS)):
   ABS[i] = 0.5 * ABS[i - 1] + gyro_bias_stability * np.random.randn(1)

##################################
####  Compute Allan Deviation  ###
##################################
output = omega + arw + ABS           # simulated sensor output
theta = np.cumsum(output) * 1 / fs   # integrate to get angles

tau, dev = allan_deviation(theta, fs=fs, out_size=200)

log_dev = np.log10(dev)
log_tau = np.log10(tau)
ddev_dtau = np.diff(log_dev) / np.diff(log_tau) 

i = np.abs(ddev_dtau).argmin()
plt.loglog(tau, dev, tau, [dev[i]] * len(tau), '--', tau[i], dev[i], 'o');

B = dev[i] * np.sqrt(2 * np.log(2) / np.pi)
print(f'Derived Input: \t{B} \nActual Input: \t{gyro_bias_stability}')

This code results in the following plot, and the derived input B. As you can see from the plot and the printed values, B is off from actual input instability by over an order of magnitude.

I'm just wondering where I went wrong in modelling the bias instability?