
import numpy as np
import matplotlib.pyplot as plt

def autocorr_formula(j0: int, i0: int, n: int) -> float:
    try:
        n = int(n)
    except:
        raise TypeError('The variable n passed into the function autocorr_formula() cannot be converted into the int type')
    if n <= 0:
        raise ValueError('The variable n passed into the function autocorr_formula() must be a positive integer')
    try:
        i = abs(int(i0))
    except:
        raise TypeError('The variable i passed into the function autocorr_formula() cannot be converted into the int type')
    try:
        j = abs(int(j0))
    except:
        raise TypeError('The variable j passed into the function autocorr_formula() cannot be converted into the int type')    
    if int(j0)*int(i0) < 0:
        raise ValueError('The integers j and i passed into the function autocorr_formula() must be of the same sign')
    elif j == i:
        return 1
    elif j > i:
        i, j = j, i
    return float((6.0 * j * np.power(n, 2) + 3.0 * np.power(n, 2) - 3.0 * n) / (np.sqrt(6.0 * i * np.power(n, 2) + 2.0 * np.power(n, 2) - 3 * n + 1) * np.sqrt(6.0 * j * np.power(n, 2) + 2.0 * np.power(n, 2) - 3 * n + 1)))

def autocov_formula(j0: int, i0: int, n: int) -> float:
    try:
        n = int(n)
    except:
        raise TypeError('The variable n passed into the function autocov_formula() cannot be converted into the int type')
    if n <= 0:
        raise ValueError('The variable n passed into the function autocov_formula() must be a positive integer')
    try:
        i = abs(int(i0))
    except:
        raise TypeError('The variable i passed into the function autocov_formula() cannot be converted into the int type')
    try:
        j = abs(int(j0))
    except:
        raise TypeError('The variable j passed into the function autocov_formula() cannot be converted into the int type')    
    if int(j0)*int(i0) < 0:
        raise ValueError('The integers j and i passed into the function autocov_formula() must be of the same sign')
    elif j == i:
        return var_formula(j=j, n=n)
    elif j > i:
        j = i
    return n * j + (n - 1.0)/2.0

def var_formula(j: int, n: int) -> float:
    try:
        n = int(n)
    except:
        raise TypeError('The variable n passed into the function var_formula() cannot be converted into the int type')
    if n <= 0:
        raise ValueError('The variable n passed into the function var_formula() must be a positive integer')
    try:
        j = abs(int(j))
    except:
        raise TypeError('The variable j passed into the function var_formula() cannot be converted into the int type')    
    return n * (j + 1.0/3.0) + 1.0/(6.0 * n) - 0.5

def autocorr_simul(j: int, i: int, data: np.ndarray, n: int) -> float:
    ro = autocorr_formula(j0=j, i0=i, n=n)
    return float(np.corrcoef(x=data[:,i], y=data[:,j])[0,1] / (1.0 - (1.0 - np.power(ro, 2))/(2.0 * data.shape[0]) + 1.0/np.power(data.shape[0],2) ))

def autocov_simul(j: int, i: int, data: np.ndarray) -> float:
    return float(np.cov(m=data[:,i], y=data[:,j], bias=False)[0,1])

def var_simul(j: int, data: np.ndarray) -> float:
    return float(np.var(a=data[:,j], ddof=1))


duration = 600
n = 6
sample = 100000
i = 29
j = 5
if i > duration // n:
    i = duration // n
if j > duration // n:
    j = duration // n
    
wiener = np.random.normal(size=((sample, (duration // n) * n - 1)))
wiener = np.hstack((np.zeros((sample,1)),wiener))
wiener = wiener.cumsum(axis=1).reshape((sample, (duration // n), n)).mean(axis=2)

print('The simulation with %d time points in total\naveraged over %d time points for the units %d and %d.' % (duration, n, j, i))
print('The time points and units are indexed from zero.')
print('\t\t\t\tby formula\tMonte Carlo')
print('Autocorrelation coefficient:\t%10.7f\t%10.7f' % (autocorr_formula(i0=i, j0=j, n=n), autocorr_simul(i=i, j=j, data=wiener, n=n)))
print('Autocovariance:\t\t\t%.7f\t%.7f' % (autocov_formula(i0=i, j0=j, n=n), autocov_simul(i=i, j=j, data=wiener)))
print('Variance of the unit %3d:\t%.7f\t%.7f' % (j, var_formula(j=j, n=n), var_simul(j=j, data=wiener)))

The simulation with 600 time points in total
averaged over 6 time points for the units 5 and 29.
The time points and units are indexed from zero.
				by formula	Monte Carlo
Autocorrelation coefficient:	 0.4368816	 0.4364241
Autocovariance:			32.5000000	32.5758148
Variance of the unit   5:	31.5277778	31.6388647

t_i , t_j	are the indices of the time points equally spaced in time
i, j	are the indices of the averaged time units
n	is the number of time points in each averaged unit
σ²	is scaled to the distance between neighboring time points

Monte Carlo verification of analytically derived autocovariance and autocorrelation matrices for an averaged Wiener process¶

Notation¶

Autocovariance matrix of the discrete-time Wiener process¶

Autocorrelation matrix of the discrete-time Wiener process¶

Autocovariance matrix of the averaged Wiener process¶

Autocorrelation matrix of the averaged Wiener process¶