Sample from a sum of two RVs by convolution of their samples via fourier transform in Excel

Question

I apologies if this is off topic for this site, but I am trying to figure out (and not having much luck) how to get an empirical distribution function of a sum of two random variables given two samples of equal length (say 256) from their respective distributions. I understand that the distribution is given by a convolutions of their distributions, which can be computed as a product of their characteristic functions.

So, to test, I take two time series of length 256 each from standard normal distribution, compute their Fourier transforms using Fourier transform tool from Data Analysis toolpack in Excel, multiply the results using complex multiplication function IMPRODUCT, and then apply inverse Fourier transform. However, the resulting time series does not have a variance of 2 (it is more like ~240) which should be the case for the sum of two standard normal RVs. It might be some normalization issue, but I am not sure how and why it comes up.

Also, I am not sure if I am supposed to be doing zero-padding of the two original time series (say +256 zeros for each). I have read that if one wants to perform regular convolutions (instead of a "circular" one, whatever that is), one have to do this padding with zeros. This padding, however, seems to make the final result also double the length, but I would like to get the estimated distribution of the sum of the same length as the inputs, so not sure what needs to be done here.

Add 1

For example, given two identical samples of size 9 from uniform distribution

$(0,0.125,0.250,0.375,0.50,0.625,0.750,0.875,1.0)$,

their convolution would produce 81 pairs with values

$(0,0.125,0.250,0.375,0.50,0.625,0.750,0.875,1.0,1.125,1.250,1.375,1.50,1.625,1.750,1.875,2.0)$

with respective counts

$(1,2,3,4,5,6,7,8,9,8,7,6,5,4,3,2,1)$

which gives the triangular distribution. What I am looking to do, is to produce this without doing these convolution sums, but by using Fourier transform.

Add 2

This is what Wikipedia says on the calculation of empirical characteristic function.

If $X_{i},i=1,2,\ldots ,n$ are i.i.d. observations, then the empirical characteristic function $\varphi _{n}(t)$ is defined as

$$\varphi _{n}(t)={\frac {1}{n}}\sum _{j=1}^{n}\exp(itX_{j}).\tag{1}$$

And this is what it says about application of the characteristic function to sums of random variables.

In particular, $\varphi _{X+Y}(t) = \varphi _{X}(t) \varphi _{X}(t)$, i.e. the characteristic function of the sum of RVs is the product of their respective characteristic functions.

So, from the above two, it seems that one can construct characteristic functions directly from the samples and multiply them together to get a characteristic function of the sums. The only bit left is to get back into the "sample space". That was my understanding.

Add 3

A discreet Fourier transform (DFT) of a sequence is

$$\hat f_k=\sum _{j=0}^{n-1} f_j \exp(i2\pi k j/n).\tag{2}$$

So, if we take the probability of observing each value in the sample on $n$ observations the same, ie $f_j = 1/n$, then the difference between empirical characteristic function and DFT is that $X_j$ in the former is replaced with $2\pi j /n$ in the latter, so the characteristic function looks like some kind of "non-uniform" DFT. And this is where I seem to get stuck, I don't know how to compute this non-uniform DFT. Also, the discreet $t$ points of the transform need to be the same for both samples as I would need to take the product of the two empirical characteristic functions at those points.

Answer 1

I still think you're completely confused about what you're trying to do, so I'm going to go back and assume your question is:

Given two data sets $x_n, n=0,\ldots, N-1$ and $y_n, n=0,\ldots, N-1$, how do I compute the probability density function (PDF) of $x_n + y_n$ using Excel?

and try to answer that.

Finding the PDF of $x_n$

Finding the empirical PDF of a data sequence is a little tricky. The most straightforward way is to find the histogram of the data and normalize it so that it sums to 1.

In pseudocode:

h_x = histogram(x)
pdf_x = h_x/sum(h_x)

It's not Excel, but this Google Sheet shows how that can be achieved.

The columns are:

• x : the random sequence.
• bin start : the start of each bin (next one is the end of this bin)
• count : the number of samples of x that fall between the bin start and bin end.
• pdf_hat : a stab at the empirical PDF using the histogram.

The same approach can be used for finding the PDF of $y$.

Finding the PDF of $x_n + y_n$

As you say, the PDF of the sum is the convolution of the two individual PDFs (provided $x_n$ and $y_n$ are independent). See section 7.2 here.

Unfortunately, Google Sheets doesn't have an FFT function, but you could do it element by element, which is what I've tried to show in the columns:

• index : index of the PDF estimates.
• Realization 1 : the pdf hat values for one realization.
• Realization 2 : the pdf hat values for another realization.
• Convolution : the element-by-element convolution.

The resulting graph is an estimate of the PDF of the sum, which approximates the expected triangle function quite well.

How to do better?

Generally, a better way to estimate the PDF of a random variable is to use kernel density estimators. Once you've got the individual PDF estimates then you can just convolve them. Or, if you have the summed data, just apply the kernel density estimator directly to that data.

Also, as you say, a better way to do convolution may be to use the FFT. To avoid "time aliasing" you need to zero-pad your signals (PDFs in this case) to be at least $N+M-1$ long (where one signal is length $N$ and the other $M$).