We have defined the covariance and correlation functions of a stationary process in Section 1.2.2. We now illustrate how to obtain the covariance and the correlation functions of a given model using this package. The functions

CovarianceFunction[model, h] and CorrelationFunction[model, h]

give, respectively, the covariance and the correlation functions of the given model up to lag h, that is, {(0), (1), ... , (h)} and {(0), (1), ... , (h)}. The code for each function solves internally a set of difference equations obtained by multiplying both sides of (2.3) by X_t-k (k=0, 1, ... ) and taking expectations. For mathematical details see Brockwell and Davis (1987), p. 92.

We begin by considering the behavior of the covariance and correlation functions of AR models.

These are the covariances at lags 0, 1, ..., 4. Note that the first entry in the output is (0), the variance of the series. To get the correlation function we need to divide the above result by its first entry (0). This can be done explicitly as follows.

We may also obtain the correlation function directly.

Both CorrelationFunction and CovarianceFunction use the function StationaryQ to check the stationarity of the model before computing the correlation or covariance function. If the model is manifestly nonstationary, the covariance or the correlation is not calculated.

When symbolic coefficients are used, StationaryQ will not give True or False and the covariance and the correlation functions are calculated assuming the model is stationary, as we have seen in Example 2.5.

Example 2.5 shows that the correlation function of an AR(1) process decays exponentially for a stationary model, and when ₁<0 it oscillates between positive and negative values. In order to visualize the "shape" of the correlation function, we plot it as a function of the lag. Since the correlation function is a list of discrete values, we use the Mathematica function ListPlot to plot it. ListLinePlot could be used instead if a line plot through the values is desired.

The output of CorrelationFunction or CovarianceFunction corresponds to lags 0, 1, 2, ..., while ListPlot assumes x coordinates of 1, 2, ... if x coordiantes are not explictly given. To match the lags with the correct correlation terms, we can either drop the first entry of the output of CorrelationFunction and plot {(1), (2), ... , (h)} using ListPlot, plot the output using ListPlot with an additional DataRange option, or input the correlation data to ListPlot in the form {{0, (0)}, {1, (1)}, ... , {h, (h)}}. This form can be obtained by using Transpose.

The output can then be used as the argument of ListPlot. It is convenient to define a function that takes the output of CorrelationFunction and does the appropriate plot so the data does not need to be manually processed every time we plot the correlation function. We will choose to use the DataRange option to ListPlot and define the function plotcorr as follows.

The first argument in plotcorr is the output of CorrelationFunction and the second is the usual set of options for ListPlot. Next we display the two typical forms of the correlation function of an AR(1) process, one for positive ₁(=0.7) and one for negative ₁(=-0.7) using the function plotcorr. (Since the correlation function of a univariate time series is independent of the noise variance, we set ²=1.)

Note that if we do not need the correlation function for other purposes we can directly include the calculation of the correlation function inside the function plotcorr as shown below.

We have given the option AxesLabel -> {"k", "(k)"} to label the axes. We can also specify other options of ListPlot. (To find out all the options of ListPlot use Options[ListPlot].) For example, if we want to join the points plotted, then the option Joined -> True should be given, and if we want to label the plot we use PlotLabel -> "label" as in the following example.

The way the correlation function decays is intimately related to the roots of (x)=0. Complex roots give rise to oscillatory behavior of the correlation function as we observe in this example. For the explicit expression of the covariance function (k) in terms of the zeros of (x), see Brockwell and Davis (1987), Section 3.3.

Next we study the behavior of MA models. Recall that the covariance function (k) of an ARMA process is calculated by multiplying both sides of (2.3) by X_t-k and computing expectations. Note that for an MA(q) process when k>q there is no overlap on the right-hand side of (2.3). Thus (k)=0 ((k)=0) for k>q. This is characteristic of the MA correlation function, and it is, in fact, often used to identify the order of an MA process, as we shall see in Section 1.5.

We see that (k)=0 for k>2.

In fact, the correlation function of an MA model can be easily worked out analytically (see Brockwell and Davis (1977), p. 93). In particular, when an MA(q) model has equal weights (i.e., ₁=₂=... =_q=), the correlation function is given by (0)=1, (k)=0 for k>q, and (k)=(+²(q-k))/(1+²q) for 0<k≤q. In particular, when =₀=1, the correlation function is (k)=(1+q-k)/(1+q) for k≤q, a straight line with slope -1/(1+q). (For convenience, we define ₀=1 so that the MA polynomial can be written as (x)=_ixⁱ. Similarly we write the AR polynomial as (x)=_ixⁱ with ₀=1.)

Note that we have avoided typing t1 eight times by using Table[t1, {8}] to generate the eight identical MA coefficients. To get the correlation function for =1 we can simply substitute t1=1 in the above expression using % /. t1 -> 1.

To emphasize the discrete nature of the correlation function some people prefer to plot the correlation function as discrete lines joining the points {i, 0} and {i, (i)} for i=0, 1, ... , h. It is easy to implement this type of plot in Mathematica, via ListPlot with a Filling option.

Using CorrelationFunction we can generate the correlation functions of different models; plotting the correlation functions enables us to develop intuition about different processes. The reader is urged to try a few examples.