# $i^{\text{th}}$-dimensional autocorrelation function

I am referring to the work of Stephen A. Billings on "Identification of a class of nonlinear systems using correlation analysis" from the year 1978, where it is mentioned that the $$i^{\text{th}}$$ dimensional autocorrelation function for a zero mean white Gaussian process $$x(t)$$ with spectral density of 1 W/cycle is as follows: $$$$\overline{x(t_1)x(t_2)\dots x(t_i)} = \begin{cases}0,& \text{odd } i\\ \sum_i{ \prod_{n\ne m} \delta(t_n-t_m)},& \text{even } i\\ \end{cases}$$$$

I would like to use this result in my work and check if this is true. How do I create a signal with prescribed spectral density in MATLAB? and how do I perform autocorrelation for $$i > 2$$ in MATLAB?

The statement is provably true:

1. $$x(t)$$ is zero mean and Gaussian with variance 1.
2. $$E[x(t_i)x(t_j)] = \begin{cases} 0 \mbox{ for } t_i \not= t_j \\ 1 \mbox{ for } t_i = t_j\end{cases}$$

For $$E[ x(t_1) x(t_2) \ldots x(t_N)]$$

1. For $$N$$ odd: There will always be a "left over" $$E[x(t_N)] = 0$$ multiplying the result, so it's always zero.

2. For $$N$$ even: There will be a non-zero value only if $$t_i = t_j$$ for distinct $$i$$ and $$j$$. For the specific case where all times are equal,

$$E[X^{2n}]=(2n-1)!!\sigma^{2n}$$

where !! represents the double factorial function (not two applications of the factorial). For each pair of $$t_i$$, $$t_j$$ that are equal, it will generate another 1 (since $$\sigma = 1$$ here).

Note that the paper uses time averages and the above uses ensemble averages. The paper assumes ergodicity of the time series, so the two are equivalent.

Much as it goes against my better judgement, here's some code that does this. The output validates the "proof":

For 2: should be 1 Calculated: 0.9946066112024137
For 3: should be 0 Calculated: 0.014338673370010435
For 4: should be 3 Calculated: 3.0088602833619134
For 5: should be 0 Calculated: 0.12289469288233784
For 6: should be 15 Calculated: 15.20500109670608
For 7: should be 0 Calculated: 1.6572818702690346


### Code Below

import numpy as np
import math
from matplotlib import pyplot as plt

N = 1000
Nruns = 100000

def do_simulation(M):
indices = 11*np.ones(M, dtype=int)
sum_product = 0
for idx in np.arange(Nruns):
x = np.random.normal(size=N)
sum_product = sum_product + math.prod(x[indices])

return sum_product/Nruns

def number_should_be(M):
if (M % 2 == 0):
return math.prod(range(M-1, 0, -2))
else:
return 0

print("For 2: should be " + str(number_should_be(2)) + " Calculated: " + str(do_simulation(2)))
print("For 3: should be " + str(number_should_be(3)) + " Calculated: " + str(do_simulation(3)))
print("For 4: should be " + str(number_should_be(4)) + " Calculated: " + str(do_simulation(4)))
print("For 5: should be " + str(number_should_be(5)) + " Calculated: " + str(do_simulation(5)))
print("For 6: should be " + str(number_should_be(6)) + " Calculated: " + str(do_simulation(6)))
print("For 7: should be " + str(number_should_be(7)) + " Calculated: " + str(do_simulation(7)))

• Thank you for the detailed explanation. Do you have an answer why the spectral density of the signal should be 1 W/cycle? And How is this achieved in your Python code? Does the random.normal function automatically returns a signal with spectral density 1? Commented Jan 16 at 12:17
• Also, I would like to ask whether you are acquainted with this paper? The paper is pretty old, some of the references are not available, and the explanation in the paper is not quite comprehensible. I have a Wiener-Hammerstein system (dynamic linear block followed by static nonlinear block and by dynamic linear block) with input $u$ and output $y$. I would like to understand, how the parameters of these individual blocks are estimated in this paper. Commented Jan 16 at 14:42
• @Neuling I'm just assuming that the 1 W/cycle means a standard deviation of 1. That's achieved by np.random.normal(size=N) which has a mean of 0 and a standard deviation of 1 as the default values. Your second comment: that sounds like a whole other question.
– Peter K.
Commented Jan 21 at 15:43