Stochastic Processes

B.4 Stochastic Processes

This appendix will provide a very brief overview of stochastic (or random) processes. If the reader finds extra time, he/she can find many good textbooks about them.

Stochastic or random processes are a powerful formalism for studying signals that incorporate some randomness. While the outcome of a random variable is a number, the outcome of a stochastic process is a time-series, each one called a realization. Hence, when using stochastic processes, it is possible to calculate “time” statistics of a single realization. Besides, one can choose a given time instant and calculate “ensemble” statistics over all realizations at that specific time. These two options may be confusing at first, but are essential. To get intuition on them, the reader is invited to observe the following simple example.

Example B.3. Simple example to explain time and ensemble averages in a stochastic process. Assume a vector [2, 0, 2, 3, 3, 2, 3] where each element represents the outcome of a random variable $X$ indicating the number of victories of a football team in four matches (one per week), along a seven-months season (for example, the team won 2 games in the first month and none in the second). One can calculate the average number of victories as $𝔼 [X] \approx 2.14$ , the standard deviation $σ \approx 1.07$ and other moments of a random variable. This vector could be modeled as a finite-duration random signal $x [n]$ , with Figure B.3 depicting its graph.

Figure B.3: Example of a finite-duration random signal.

Assume now there are 5 vectors, one for each year. A discrete-time random process can be used to model the data, with each vector (a time series) corresponding to a realization of the random process. For example, following the Matlab convention of organizing different time series in columns (instead of rows), the tabular data below represents the five random signals:

victories = [
     2     3     0     0     1
     0     1     3     2     1
     2     1     4     0     2
     3     4     3     1     2
     3     3     4     3     4
     2     1     1     2     3
     3     4     4     2     4]

Each random signal (column) can be modeled as a realization of a random process. Figure B.4 shows a graphical representation of the data with the five realizations of the example. One should keep in mind that a random process is a model for generating an infinite number of realizations but only 5 are used here for simplicity.

Figure B.4: Example of five realizations of a discrete-time random process.

A discrete-time random process will be denoted by $X [n]$ and one of its realization as $x [n]$ . For the sake of explaining this example, let us define the notation $x_{i} [n]$ to indicate the $i$ -th realization of $X [n]$ . Because $x_{i} [n]$ is a random signal, one can calculate for the given example, $𝔼 [x_{1} [n]] \approx 2.14$ and other statistics taken over time. But a random process is a richer model than a random signal and other statistics can be calculated or, equivalently, more questions can be asked.

Besides asking for statistics of one realization taken over “time”, it is important to understand that when fixing a given time instant $n_{0}$ , the values that the realizations assume are also outcomes of a random variable that can be denoted by $X [n_{0}]$ . In the example, assuming $n_{0} = 4$ (i.e., assume the fourth month), one can estimate the average $𝔼 [X [4]] = 2.6$ . Similarly, $𝔼 [X [6]] = 1.8$ . Figure B.5 graphically illustrates the operation of evaluating a random process at these two time instants. Differently than the previous time averages, these “ensemble” statistics are obtained across realizations with fixed time instants.

Figure B.5: Example of evaluating a random process at time instants $n = 4$ and $n = 6$ , which correspond to the values of two random variables $X [4]$ and $X [6]$ . There are only five pairs and each occurs once, hence a 0.2 estimated probability for all points.

The continuous-time random process $X (t)$ is similar to its discrete-time counterpart. Taking two different time values $t$ and $s$ leads to a pair of random variables, which can then be fully (statistically) characterized by a joint pdf function $f (X (t), X (s))$ . This can be extended to three or more time instants (random variables). Back to the discrete-time example, one can use a normalized two-dimensional histogram to estimate the joint PMF of $X [4]$ and $X [6]$ as illustrated in Figure B.6.

Figure B.6: Example of a joint pdf of the continuous random variables $X [4]$ and $X [6]$ .

This example aimed at distinguishing time and ensemble averages, which is fundamental to understand stochastic processes. $□$

Correlation Function and Matrix

Second-order moments of a random process involve statistics taken at two distinct time instants. Two of them are the correlation and covariance functions, which are very useful in practice.

As discussed in Section 1.10, the correlation function is an extension of the correlation concept to random signals generated by stochastic processes. The autocorrelation function (ACF)

R_{X} (s, t) = 𝔼 [X (s) X (t)]

(B.31)

is the correlation between random variables $X (s)$ and $X (t)$ at two different points $s$ and $t$ in time of the same random process. Its purpose is to determine the strength of relationship between the values of the signal occurring at two different time instants. For simplicity, it is called $R_{X} (s, t)$ instead of $R_{XX} (s, t)$ .

A discrete-time version of Eq. (B.31) is

R_{X} [n_{1}, n_{2}] = 𝔼 [X [n_{1}] X [n_{2}]],

(B.32)

where $n_{1}$ and $n_{2}$ are two “time” instants.

Specially for discrete-time processes, it is useful to organize autocorrelation values in the form of an autocorrelation matrix.

Example B.4. Calculating and visualizing the autocorrelation matrix (not assuming stationarity). Stationarity is a concept that will be discussed in the sequel (Appendix B.4.0.0). Calculating and visualizing autocorrelation is simpler when the process is stationary. But in this example, this property is not assumed and the autocorrelation is obtained for a general process.

Listing B.1 illustrates the estimation of $R_{X} [n_{1}, n_{2}]$ for a general discrete-time random process.³

Listing B.1: MatlabOctaveFunctions/ak_correlationEnsemble.m

1function [Rxx,n1,n2]=ak_correlationEnsemble(X, maxTime) 
2% function [Rxx,n1,n2]=ak_correlationEnsemble(X, maxTime) 
3%Estimate correlation for non-stationary random processes using 
4%ensemble averages (does not assume ergodicity). 
5%Inputs:  - X, matrix representing discrete-time process realizations 
6%           Each column is a realization of X. 
7%         - maxTime, maximum time instant starting from n=1, which 
8%           corresponds to first row X(1,:). Default: number of rows. 
9%Outputs: - Rxx is maxTime x maxtime and Rxx[n1,n2]=E[X[n1] X*[n2]] 
10%         - n1,n2, lags to properly plot the axes. The values of 
11%           n1,n2 are varied from 1 to maxTime. 
12[numSamples,numRealizations] = size(X); %get number of rows and col. 
13if nargin==1 
14    maxTime = numSamples; %assume maximum possible value 
15end 
16Rxx = zeros(maxTime,maxTime); %pre-allocate space 
17for n1=1:maxTime 
18    for n2=1:maxTime %calculate ensemble averages 
19        Rxx(n1,n2)=mean(X(n1,:).*conj(X(n2,:))); %ensemble statistics 
20        %that depends on numRealizations stored in X 
21    end 
22end

Figure B.7 was created with the command ak_correlationEnsemble(victories) and shows the estimated correlation for the data in Figure B.4 of Example B.3. The data cursor indicates that $R_{X} [1, 1] = 2.8$ because for $n = 1$ the random variables of the five realizations are $[2, 3, 0, 0, 1]$ . In this case, $R_{X} [1, 1] = 𝔼 [X {[1]}^{2}] \approx (4 + 9 + 1) ∕ 5 = 2.8$ . Similar calculation leads to $R_{X} [1, 7] = 𝔼 [X [1] X [7]] \approx 4.4$ because in this case the products victories(1,:).*victories(7,:) are [6, 12, 0, 0, 4], which have an average of 4.4.

Figure B.7: Correlation for data in matrix victories.

Figure B.8 presents an alternative view of Figure B.7, where the values of the z-axis are indicated in color. The two datatips indicate that $R_{x} [1, 2] = 0.8$ and $R_{x} [3, 5] = 6.6$ (these are the values identified by Index for the adopted colormap, and the associated color is also shown via its RGB triple).

Figure B.8: Version of Figure B.7 using an image.

For real-valued processes, the autocorrelation is symmetric, as indicated in Figure B.7 and Figure B.8, while it would be Hermitian for complex-valued processes. $□$

Two time instants: both absolute or one absolute and a relative (lag)

When characterizing second-order moments of a random process, instead of using two absolute time instants, it is sometimes useful to make one instant as a relative “lag” with respect to the other. For example, take two absolute instants $t = 4$ and $s = 7$ in Eq. (B.31). Alternatively, the same pair of instants can be denoted by considering $t = 4$ as the absolute time that provides the reference instant and a lag $τ = s - t = 3$ indicating that the second instant is separated by 3 from the reference.

Using this scheme, Eq. (B.31) can be written as

R_{X} (t, τ) = 𝔼 [X (t + τ) X (t)] .

(B.33)

Similarly, Eq. (B.32) can be denoted as

R_{X} [n, l] = 𝔼 [X [n] X [n + l]],

(B.34)

with $n_{1} = n$ and $n_{2} = n + l$ , where $l$ is the lag in discrete-time.

Example B.5. Contrasting two different representations of autocorrelations. Figure B.7 is a graph that corresponds to Eq. (B.32), with two absolute time instants. When using the alternative representation of Eq. (B.34), with a lag, the resulting matrix is obtained by rearranging the elements of the autocorrelation matrix. Figure B.9 compares the two cases. Note that the name autocorrelation matrix is reserved for the one derived from Eq. (B.32) (and Eq. (B.31)) that is, for example, symmetric.

Figure B.9: Comparison between the 3-d representation of the autocorrelation matrix in Figure B.7 and the one using lags as in Eq. (B.34).

Figure B.10 is an alternative view of Figure B.9, which makes visualization easier. It better shows that Eq. (B.34) leads to shifting the column of the autocorrelation matrix (Figure B.7) corresponding to each value $n = n_{1}$ of interest.

Figure B.10: Comparison between autocorrelation representations using images instead of 3-d graphs as in Figure B.9.

Properties such as the symmetry of the autocorrelation matrix may not be so evident from the right-side representation of Figure B.10 as in the left-side representation, but both are useful in distinct situations. $□$

Covariance matrix

Similar to autocorrelation functions, it is possible to define covariance functions that expand the definition for random variables such as Eq. (1.48) and take as arguments two time instants of a random process. For example, for a complex-valued discrete-time process the covariance is

c_{X} [n, l] = 𝔼 [(X [n] - μ [n]) {(X [n + l] - μ [n + l])}^{*}] .

(B.35)

Stationarity, cyclostationarity and ergodicity

A random process is a powerful model. Hence, an important question is what information is necessary to fully characterize it. Take the $7 \times 5$ matrix in Example B.3: this data does not fully characterize the random process because it consists of only five realizations (tossing a coin five times does not allow to estimate its probability of heads robustly). We would need an infinite number of realizations if we insist in characterizing it in a tabular form. The most convenient alternative is to describe it through pdfs. We should be able to indicate, for all time $t$ , the pdf of $X (t)$ . Besides, for each time pair $(t, s)$ , the joint pdf $f (X (t), X (s))$ should be known. The same for the pdf $f (X (t), X (s), X (r))$ of each triple $(t, s, r)$ and so on. It can be noted that a general random process requires a considerable amount of statistics to be fully described. So, it is natural to define particular and simpler cases of random processes as discussed in the sequel.

An important property of a random process is the stationarity of a statistics, which means this statistics is time-invariant even though the process is random. More specifically, a process is called $n$ -th order stationary if the joint distribution of any set of $n$ of its random variables is independent of absolute time values, depending only on relative time.

A first-order stationary process has the same pdf $f (X (t)), \forall t$ (similar is valid for discrete-time). Consequently, it has a constant mean

𝔼 [X (t)] = μ, \forall t,

(B.36)

constant variance and other moments such as $𝔼 [X^{4} (t)]$ .

A second-order stationary process has the same joint pdf $f (X (t), X (s)) = f (X (t + τ), X (s + τ)), \forall t, s$ . Consequently, its autocorrelation is

R_{X} (τ) = 𝔼 [X (t + τ) X (t)],

(B.37)

where the time difference $τ = s - t$ is called the lag and has the same unit as $t$ in $X (t)$ . Note that, in general, $R_{X} (t, τ)$ of Eq. (B.33) depends on two parameters: $t$ and $τ$ . However, for a stationary process, Eq. (B.37) depends only on the lag $τ$ because, given $τ$ , $R_{X} (t, τ) = R_{X} (τ)$ is the same for all values of $t$ .

Similarly, if a discrete-time process is second-order stationary, $R_{X} [n, l]$ Eq. (B.34) can be simplified to

R_{X} [l] = 𝔼 [X [n] X [n + l]],

(B.38)

where $l = n_{2} - n_{1}$ is the lag.

A third-order stationary process has joint pdfs $f (X (t_{1}), X (t_{2}), X (t_{3}))$ that do not depend on absolute time values $(t_{1}, t_{2}, t_{3})$ and so on. A random process $X (t)$ is said to be strict-sense stationary (SSS) or simply stationary if any joint pdf $f (X (t_{1}), \dots, X (t_{n}))$ is invariant to a translation by a delay $τ$ . A random process with realizations that consist of i. i. d. random variables (called an i. i. d. process) is SSS. A random process that is not SSS is referred to as a non-stationary process even if some of its statistics have the stationarity property.

The SSS process has very stringent requirements, which are hard to verify. Hence, many real-world signals are modeled as having a weaker form of stationarity called wide-sense stationarity.

Basically, a wide-sense stationary (WSS) process has a mean that does not vary over time and an autocorrelation that does not depend on absolute time, as indicated in Eq. (B.36) and Eq. (B.37), respectively. Broadly speaking, wide-sense theory deals with moments (mean, variance, etc.) while strict-sense deals with probability distributions. Note that SSS implies WSS but the converse is not true in general, with Gaussian processes being a famous exception.

For example, the process corresponding to Figure B.7 is not WSS: its autocorrelation does not depend only on the time difference $τ$ .

Example B.6. The autocorrelation matrix of a WSS process is Toeplitz. For a WSS process, elements $R_{X} (i, j)$ of its autocorrelation matrix depend only on the absolute difference $| i - j |$ and, therefore, the matrix is Toeplitz. In this example, the function toeplitz.m is used to visualize the corresponding autocorrelation matrix as an image. For example, the command Rx=toeplitz(1:4) generates the following real-valued matrix:

 Rx=[1     2     3     4
     2     1     2     3
     3     2     1     2
     4     3     2     1].

Usually, the elements of the main diagonal correspond to $R_{X} (0) = 𝔼 [X^{2} [n]]$ .

To visualize a larger matrix, Figure B.11 was generated with Rxx=toeplitz(14:-1:4).

Figure B.11: Representation of a WSS autocorrelation matrix that depends only on the lag $l$ .

Given the redundancy in Figure B.11, one can observe that for a WSS, it suffices to describe $R_{X} (l)$ for each value of the lag $l$ . $□$

A cyclostationary process is a non-stationary process having statistical properties that vary periodically (cyclically) with time. In some sense it is a weaker manifestation of a stationarity property.

A given statistics (variance, autocorrelation, etc.) of a random process is ergodic if its time average is equal to the ensemble average. For example, an i. i. d. random process is ergodic in the mean.⁴ And, loosely speaking, a WSS is ergodic in both mean and autocorrelation if its autocovariance decays to zero.⁵

A process could be called “ergodic” if all ensemble and time averages are interchangeable. However, it is more pedagogical to consider ergodicity as a property of specific statistics, and always inform them.⁶ Even non-stationary or cyclostationary processes can be ergodic in specific statistics such as the mean. Figure B.12 depicts the suggested taxonomy of random processes using a Venn diagram.

Figure B.12: Suggested taxonomy of random processes. SSS are sometimes called stationary and all processes that are not SSS are called non-stationary.

Example B.7. Polar and unipolar bitstreams modeled as WSS processes. This example is based on processes used in digital communications to represent sequences of bits.

Figure B.13 illustrates two processes for which the main diagonal $R_{X} [1, 1] = R_{X} [2, 2] = R_{X} [3, 3] \dots$ has the same values. The process (a) (left-most) is called a unipolar line code and consists of sequences with equiprobable values 0 and 1, while for the second process (called polar) the values are $- 1$ and 1.

For the unipolar and polar cases, the theoretical correlations $R_{X} [n_{1}, n_{2}]$ for $n_{1} = n_{2}$ are 0.5 and 1, respectively. The result is obtained by observing that $R_{X} [n_{1}, n_{2}] = R_{X} [n_{1}, n_{1}] = 𝔼 [X {[n_{1}]}^{2}]$ and for the unipolar case, $X {[n_{1}]}^{2}$ is 0 or 1. Because the two values are equiprobable, $𝔼 [X {[n_{1}]}^{2}] = 0.5$ . For the polar case, $X {[n_{1}]}^{2} = 1$ always and $𝔼 [X {[n_{1}]}^{2}] = 1$ . The values of $R_{X} [n_{1}, n_{2}]$ for $n_{1} \neq n_{2}$ are obtained by observing the values of all possible products $X [n_{1}] X [n_{2}]$ . For the unipolar case, the possibilities are $0 \times 0 = 0$ , $0 \times 1 = 0$ , $1 \times 0 = 0$ and $1 \times 1 = 1$ , all with probability 1/4 each. Hence, for this unipolar example, $𝔼 [X [n_{1}] X [n_{2}]] = (3 \times 0 + 1 \times 1) ∕ 4 = 0.25, n 1 \neq n 2$ . For the polar case, the possibilities are $- 1 \times - 1 = 1$ , $- 1 \times 1 = 1$ , $1 \times - 1 = - 1$ and $1 \times 1 = 1$ , all with probability 1/4 each. Hence, for this polar example, $𝔼 [X [n_{1}] X [n_{2}]] = (2 \times - 1 + 2 \times 1) ∕ 4 = 0, n 1 \neq n 2$ .

In summary, for the polar code $R_{X} [n_{1}, n_{2}] = 1$ for $n_{1} = n_{2}$ and 0 otherwise. For the unipolar code, $R_{X} [n_{1}, n_{2}] = 0.5$ for $n_{1} = n_{2}$ and 0.25 otherwise.

Figure B.13: Correlation for random sequences with two equiprobable values: a) values are 0 and 1 (unipolar), b) values are $- 1$ and 1 (polar).

A careful observation of Figure B.13 indicates that two dimensions are not necessary and it suffices to describe $R_{X} [l]$ as in Eq. (B.34). As discussed previously, in the case of WSS processes, the correlation depends only on the difference between the “time” instants $n_{1}$ and $n_{2}$ . Figure B.14 emphasizes this aspect using another representation for the polar case in Figure B.13.

Figure B.14: Alternative representation of the correlation values for the polar case of Figure B.13.

Because the polar and unipolar processes that were used to derive Figure B.13 are WSS, Listing B.2 can be used to take advantage of having $R_{X} [n_{1}, n_{2}]$ depending only on the lag and convert this matrix to an array $R_{X} [l]$ .

Listing B.2: MatlabOctaveFunctions/ak_convertToWSSCorrelation.m

1function [Rxx_tau, lags] = ak_convertToWSSCorrelation(Rxx) 
2% function [Rxx_tau,lags]=ak_convertToWSSCorrelation(Rxx) 
3%Convert Rxx[n1,n2] into Rxx[k], assuming process is 
4%wide-sense stationary (WSS). 
5%Input: Rxx -> matrix with Rxx[n1,n2] 
6%Outputs: Rxx_tau[lag] and lag=n2-n1. 
7[M,N]=size(Rxx); %M is assumed to be equal to N 
8lags=-N:N; %lags 
9numLags = length(lags); %number of lags 
10Rxx_tau = zeros(1,numLags); %pre-allocate space 
11for k=1:numLags; 
12    lag = lags(k); counter = 0; %initialize 
13    for n1=1:N 
14        n2=n1+lag; %from: lag = n2 - n1 
15        if n2<1 || n2>N %check 
16            continue; %skip if out of range 
17        end 
18        Rxx_tau(k) = Rxx_tau(k) + Rxx(n1,n2); %update 
19        counter = counter + 1; %update 
20    end 
21    if counter ~= 0 
22        Rxx_tau(k) = Rxx_tau(k)/counter; %normalize 
23    end 
24end

Figure B.15 was generated using Listing B.2 for the data in Figure B.13.

Figure B.15: One-dimensional ACF $R_{X} [l]$ for the data corresponding to Figure B.13 (unipolar and polar codes).

If ergodicity can be assumed, a single realization of the process can be used to estimate the ACF. Figure B.16 was generated using the code below, which generates a single realization of each process:

1numSamples=1000; 
2Xunipolar=floor(2*rand(1,numSamples)); %r.v. 0 or 1 
3[Rxx,lags]=xcorr(Xunipolar,9,'unbiased');%Rxx for unipolar 
4subplot(121); stem(lags,Rxx) title('Unipolar'); 
5Xpolar=2*[floor(2*rand(1,numSamples))-0.5]; %r.v. -1 or 1 
6[Rxx,lags]=xcorr(Xpolar,9,'unbiased'); %Rxx for polar code 
7subplot(122); stem(lags,Rxx) title('Polar');

Figure B.16: ACF estimated using ergodicity and waveforms with 1,000 samples for each process. The theoretical values are 0.25 and 0.5 (at $l = 0$ ) for unipolar, while for polar these values are 0 and 1.

Increasing the number of samples to numSamples=1e6 achieves estimation errors smaller than in Figure B.16. $□$

³ Note that xcorr.m in Matlab/Octave assumes the process is wide sense stationary, while ak_correlationEnsemble.m does not.

⁴ See, e. g., [Kay06], pg. 564.

⁵ See [ urlBMerm] for the requirement of a WSS to be ergodic in the mean and/or [Kay06], pg. 577.

⁶ See [ urlBMerg] for an example of a SSS process that is not ergodic. The discussion also illustrates how it is confusing when the statistics that is ergodic is not specified and the process is called “ergodic”.