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