import numpy as np, math
rng=np.random.default_rng(1414)
def samples(rho,b,trials=8000):
 e=rng.normal(size=(trials,b)); z=np.zeros_like(e); z[:,0]=e[:,0]; q=math.sqrt(1-rho*rho)
 for i in range(1,b): z[:,i]=rho*z[:,i-1]+q*e[:,i]
 return z.mean(1)-z[:,:b//2].mean(1)
print('coupled correction variance, rho=.8, b0=32')
for l in range(5):
 b=32*2**l; v=samples(.8,b).var(); print(l,b,v,v*b)
print('sample mean variance ratios vs iid, predicted long-run multiplier')
for rho in [0,.5,.8,.95]:
 vals=[]
 for b in [32,128,512]:
  e=rng.normal(size=(4000,b)); z=np.zeros_like(e); z[:,0]=e[:,0]; q=math.sqrt(1-rho*rho)
  for i in range(1,b): z[:,i]=rho*z[:,i-1]+q*e[:,i]
  vals.append(z.mean(1).var())
 print(rho, vals, 'ratio', [v/vals[0] for v in vals], 'pred', (1+rho)/(1-rho))