Exploring Data Graphically in Mathematica

Robert A. Stine
Statistics Department
The Wharton School
University of Pennsylvania
Philadelphia, PA  19104-6340

www-stat.wharton.upenn.edu/~stine

September 17, 1995
Revised November 4, 2003

Introduction

In 1995, I produced a Mathematica package that I called DataTools.  This package implemented a variety of statistical functions that were then either hard to do or not available in the standard distribution of Mathematica. Since then, Mathematica has grown, and the supplied packages have now filled the void that used to exist.

This notebook shows how to do the same things that I once implemented "from scratch" in DataTools, but uses features supplied in Version 5 of Mathematica.  The problems (and solutions) are the same, but rather than using  programs that I developed, I'm going to show how to make the most of now-included functions.  For example, Mathematica now has a nice histogram function.  It's also got boxplots.  I'll use Mathematica's version of the first, but I still like my own boxplot.  I'll show you both and let you choose.

To present these functions in a more realistic (and hopefully interesting) context, I introduce them through a small exploratory analysis of the daily returns (including dividends) on General Motors common stock during 1987-1988.  This period includes the big drop in October, 1987.  The data are in the file gm8788.dat which is available from my Web page.  The original source for this data is CRSP, the Center for Research in Security Prices.

The analysis of these data progresses in a fashion intended to introduce the functions of originally presented in the DataTools package, and so might seem a bit contrived. After reading the data, the analysis begins with simple descriptive statistics that use binning to summarize the data a month at a time. The next two sections use univariate and bivariate plots, including histograms, boxplots, and kernel density estimates.  The last 'technical' section explores how the volatility in these data tracks over time using regression.  Throughout, the emphasis is upon plots that reveal interesting features.

To use the functions that I have written, you'll need to load these functions.  They're defined in the file DataTools.m.

\(\(<< "/Users/bob/Documents/Mathematica/DataTools/DataToolsCode.m";\)\)

Reading and Writing Data

ReadData

The format of a data file seems a personal choice, given the variety of formats out there.  My function for reading data simply assumes that the data has the following layout.  ReadData expects the data file to possess a very simple format which is often used by commercial statistics packages for interchanging data.  The first line of the file is a quoted string used to document the file;  this allows you to insert a comment to remind you of the file's contents.  The second line of the file is a list of variable names separated by one or more spaces.  The remaining lines of the file are the observation records, with as many values per line as there are variable names.  For example, a file with two variables and 3 observations would appear as follows:

"A simple file example"
  Var1  Var2
  2  3
  12 na
  1    3

Missing data can be denoted by whatever symbol you choose, such as the na symbol used here or any symbolic tag such as the built-in symbol Indeterminant.  ReadData prints the documentation string.  It then defines the symbols, here Var1 and Var2.  After running ReadData on this file, you'll have two new variables, each a list with 3 values.  ReadData avoids over-writing a pre-existing symbol, as explained further below.

The data file gm8788.dat contains 5 columns of daily information: the day  month, year, price and return on General Motors common stock for the 506 trading days during 1987 and 1988.  Each line of the file (after the first two descriptive lines) contains five values that define these variables for each day.  To distinguish the names of variables read from the file from other symbols, I used upper-case characters in the name for each (for example, DAY for the day of the month).   Here are the first few lines of the data file gm8788.dat, along with the output of the ReadData function.  (You'll need to tell Mathematica where to look for the data file on your system.)

"Returns on General Motors Stock during 1987-1988"
DAY MONTH YEAR PRICE RETURN
02 1 1987 66.875 0.01325758
05 1 1987 68.5 0.02429906
06 1 1987 68.75 0.003649635
07 1 1987 69.125 0.005454545

ReadData["~/Documents/Mathematica/DataTools/gm8788.dat"]

After running this function, the named columns are defined.  A list named by the associated symbol extracted from the second line of the input file holds the data for each variable.  For example, an abbreviated list of the returns  is

Short[RETURN]

This output shows the first three daily returns, then indicates 500 values are hidden, and then shows the last three daily returns.

ReadData does not overwrite  previously defined symbols.  If some of the symbols which are to denote variables (in this example, the symbols DAY, MONTH, YEAR, PRICE and RETURN) are already defined in the active session, a warning message is printed.  This feature avoids accidentally replacing an important piece of data by new data read from the external file. For example, a second attempt to read the same file generates a message but does not define any variables.

ReadData["~/Desktop/gm8788.dat"]

To re-read data (if you've made a change to the data, for example)  use Clear to remove the old definitions of the conflicting symbols (as in Clear[DAY, MONTH, YEAR, PRICE, RETURN]) and repeat the  read.

WriteData

The function WriteData creates a file with this format by combining several column lists of values.  With the data defined, we can now turn to some basic plots and summary statistics.

WriteData[{RETURN}, "~/Desktop/test.dat",  "One column"]

Descriptive Statistics, Sequence Plots, and Binning

Descriptive Statistics

Mathematica now includes various descriptive statistics such as the mean and standard deviation. With the improvements in the underlying numerical algorithms, these run quite fast. Once upon a time (check out the original description of this package in Varian) it mattered a great deal how one defined numerical calculations in Mathematica.  With the improvements in speed of the host systems and the kernel of Mathematica itself, these concerns are pretty irrelevant now.  That said, I suggest using the built-in mean function rather than rolling your own.  Its a lot faster.

Here are some basic examples. The average and standard deviation of the returns over these two years are:

ave = Mean[RETURN]

sd = StandardDeviation[RETURN]

Other summaries include the variance, the median, and the interquartile range. They've even included the MAD, or median absolute deviation about the median.

MedianDeviation[RETURN]

Missing Values

These summaries are pretty simple unless you have missing data. For example, if we just add a symbol to the end of the returns, the averaging function no longer works.  Well, not in the sense that you might want for analyzing data.  You end up with an expression that includes the missing symbol.  Fine for symbolic algebra (it makes it fun to draw the influence curve), but not what you want if the data is simply missing.

Mean[Append[RETURN, miss]]

So, to handle missing data you need to first filter the data to remove terms that are not numbers.  The function DropNonNumeric from Statistics`DataManipulation` handles this conveniently.

Mean[DropNonNumeric[Append[RETURN, miss]]]

Care About Speed?

Before moving on, these simple functions offer some insight into the speed Mathematica.  The built in functions use the developers considerable knowledge of how to make the code run faster. Consider the time that it takes to compute the average of 100,000 uniform random numbers.

x = Table[Random[], {100000}];

The obvious way to compute the average is to use the function Plus.  Repeating the calculation 100 times takes about 2.2 seconds. (The command Plus@@x is read as the summation  x.)

First[Timing[Table[Plus @@ x, {100}]]]

In comparison, the built-in function is about 50 times faster.  You can see how they get this to happen by looking at the code that implements the various descriptive statistics.  It pays to recognize the data is a packed array.

First[Timing[Table[Mean[x], {100}]]]

Sequence Plots and Dates

A careful analysis of any time series begins with sequence plots.  The built-in function ListPlot requires some simple but tedious manipulations which my wrapper ScatterPlot provides automatically when given a single list as its argument. ScatterPlot is basically a call to ListPlot, but Scatterplot adds some touches that are important in data analysis, including a transposed argument list and setting the option PlotRange set to show all of the data.  The options of Scatterplot, such as that to add a label to the plot, are a superset of those of ListPlot.

ScatterPlot[PRICE, 
        PlotLabel -> "Daily GM Stock Price"];

In contrast to the drift in  daily stock prices, the associated returns (that is, the change in price from the previous trading day divided by the price on the previous day) show no evident trend.

ScatterPlot[RETURN, 
        PlotLabel -> "Daily GM Stock Price"];

With a little more effort, we can give the horizontal  axis a more useful set of labels that identify the time of the big drop in the prices, and the large swing in returns (as if you did not know when these happened).  The first step is to associate dates with the axis labels. Since October 19, 1987 (day number 202), is the day the stock market fell substantially, it's useful to label this specific date.

Position[RETURN, Min[RETURN]]

Having a function that forms the labels makes the calculations more easy to follow.  Here is a simple version that builds a date string.  The result is a (numerical value, date string) pair.

\(calDate[index_] := 
      With[{m = {"Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"}}, {index, 
          ToString[
              DAY[[index]]] <> " " <> 
            m[[MONTH[[
                index]]]] <> " " <> 
            ToString[
              YEAR[[
                index]]]}]\)

calDate[1]

These labels are compatible with the standard Ticks option used in plots in Mathematica. It's probably also a good idea with these extended tick marks to move the axis beneath the data and make the plot relatively longer. Also, since such date labels will crowd on the axis, the Tick option below provides just three dates. The ticks for the y-axis are determined automatically.

\(\(sp = 
        ScatterPlot[RETURN, PlotLabel -> "Daily GM Stock Return", 
          Ticks -> {calDate /@ {1, 202, 506}, 
              Automatic}, 
AxesOrigin -> {0, \(-0.25\)}, 
AspectRatio ->  .5];\)\)

We can further refine of this plot by adding, for example, grid lines based on the average and standard deviation of the returns.

m  =  Mean[RETURN];  sd  =  StandardDeviation[RETURN]; 
    Show[sp, GridLines -> {None, {m + 2 sd, m - 2 sd}}];

It appears that the outlying values have inflated the standard deviation so that more than the usual 95% of the returns fall within two standard deviations of the mean.  Using the following boolean function together with Count, we find that the actual percentage within this interval is closer to 97%.

withinQ[x_] := 
      m - 2 sd ≤ x ≤ m + 2 sd

N[Count[RETURN, x_?withinQ]/Length[RETURN]]

Binning

Simple descriptive statistics like Average and StandardDeviation become more interesting and powerful when combined with the grouping function Bin.  As with built-in functions and those that I have defined in DataTools, each of my functions has an accompanying description accessed by the query function ``?''.

?Bin

Here is a simple example of how to use Bin. In this example, each of the values in the first list that share a common value (1 or 2) are grouped.

Bin[{1, 1, 1, 2, 2, 2,  1, 2}, {1, 2, 3, 4, 5, 6,  7,  
        8}]

Binning makes it easy to summarize the daily stock returns by calendar month.  In the following examples, the trailing semi-colon  avoids a lengthly print-out on the screen.  Also, Bin will print a warning message if, as illustrated in the next command, the arguments appear in the wrong order and thereby attempt to create a substantial number of bins (the cutoff in the code is 90% of the length of the input list).

byMonth = Bin[RETURN, MONTH];

The next command groups the data by month into a list with 24 two-term lists.

monthCode  =  
        N[YEAR + ((MONTH - 1)/12)];
byMonth = Bin[monthCode, RETURN];
Length[byMonth]

The first member of each paired list in byMonth is the numerical date, and the second member is the list of associated values.  Here is the paired list for January, 1987.

First[byMonth]

With the returns grouped into monthly subsets, function mapping gives summary statistics for each month.  In the next command, I apply Mean and StandardDeviation to the returns for each month (second member of each list).  Each date (the first member)  is converted into a string, dropping the redundant characters ``19''.  Each item in the resulting summary list is a triple giving the mean, standard deviation, and month label.

summary = (
{Mean[#1[[2]]], 
StandardDeviation[#1[[2]]], 
StringDrop[
                  ToString[#1[[1]]], 
                  2]} &) /@ byMonth;
First[summary]

The format of this summary as a list triples is ideal for the plotting functions that come with Mathematica.  These functions expect data to be organized as a list of observations, with each element of the list denoting the values for one observation.  The function LabeledListPlot from Graphics`Graphics` gives a nice plot; DataTools automatically searches this package.  The default scaling, however, conceals the large outlier from October, 1987; the option PlotRange->All forces the plot to show all of the data.

LabeledListPlot[summary, 
        AxesLabel -> {"Ave Return", "SD Return"}, 
        PlotRange -> All];

To see the interesting pattern in the other 23 months, we can manually set the range.

LabeledListPlot[summary, 
        AxesLabel -> {"Ave Return", "SD Return"}, 
        PlotRange -> {{- .005,  .010}, { .005,  .025}}];

As financial theory suggests, months with higher returns also tend to have higher standard deviation.  That is, risk and return are positively correlated if we set aside October, 1987.

To compute the correlation explicitly, first transpose the list of means and standard deviations so that the columns are easily accessed.  With October, 1987, included in the data, the correlation is negative.  (The function Correlation is part of the package Statistics`MultiDescriptiveStatistics`; DataTools imports this one as well.)  Because the correlation function expects a pair of lists as arguments, not a list of pairs, we have to do some rearranging.

The first mapped functioin pulls off the first two summary elements (which are the mean and the standard deviation.  The @@ then makes the two columns of these when transposed the inputs to the correlation function.  (Yuk.  Too bad that the Correlation function is not like the plot functions and works with a list of pairs instead.  Or at least have them all work the same way.)

Correlation @@ 
      Transpose[{#[[1]], #[[2]]} &  /@  summary]

We can set aside October, 1987, (the tenth month) without losing its value by simply tagging it with a symbolic value.  The tagging is done by multiplying this observation in place using the ``C-like'' expression.

summary[[10, 1]] *= na

With this outlier set aside (its still there, just not used in this calculation), the correlation is virtually the same as before, but positive.  DropNonNumeric works on lists of pairs.

Correlation @@ 
      Transpose[
        DropNonNumeric[{#[[1]], #[[2]]} &  /@  
            summary]]

My function Bin also supports some options which are useful when the binning factor is not categorical.  These options require that the data be sorted so that the leading argument, which defines the binning factor, is in the desired order.

Options[Bin]

The default bins in this case hold either a fixed proportion of the elements of the second argument as ordered by the first argument, or group the data over intervals of a fixed width.  The leading element of the resulting bin lists is the average of the associated values defined by the first argument.  For example, the next command splits the data into two bins.

Bin[{10, 20, 30, 40, 50, 60}, {1, 2, 3, 4, 5, 6}, 
      BinPercentage -> 0.5]

To illustrate the use of this type of binning,  let's first build a numerical date by combining the day, month, and year information.  Ignoring the lengths of the months, we have

DATE = DAY + MONTH 30 + (YEAR - 1987) 360;
Short[DATE]

Below is a sequence plot of the average returns for 30 day segments. The argument to ListPlot is formed by function mapping. The first element is just the first element of each list produced by binning (the ``average date''), and the second is the average return.  As with other built-in plotting commands, we need to set the option PlotRange -> All to avoid hiding outlying portions of the data.  (You can use SetOption to change the default for the entire session.)

ListPlot[({#1[[1]], 
                Mean[#1[[2]]]} &) /@ 
          Bin[DATE, RETURN, BinWidth -> 30], 
        PlotRange -> 
          All, 
PlotLabel -> "Mean returns by \
month"];

Distribution Plots

Histograms

The Histogram function from Graphics`Graphics` draws this familiar summary of the distribution of a set of values.  Continuing our analysis of the GM returns, the default bin width is small because of the extreme outliers noted in the initial time series plots.  The option HistogramCategories controls the number of bins.  The option HistogramScale controls the area of the plot; setting this to 1 gives a density.

Histogram[RETURN, HistogramCategories -> 30,  
        HistogramScale -> 1];

The built-in function Select makes it easy to filter out the extremes, which in this example are associated with the large volatility of the stock market in October, 1987.  The in-line function identifies the returns that fall in this range.  (The curious labeling of the x-axis shows that recognizing zero is not so easy as it seems.)

\(\(hist = 
        Histogram[Select[RETURN, Abs[#1] < 0.1 &], 
          PlotLabel -> "Filtered GM Return", 
HistogramCategories -> 30, \ 
          HistogramScale -> 1];\)\)

Kernel Density Estimates

The histogram conveys an inadequate impression of the distribution of the underlying population.  Kernel density estimates (Silverman 1986) eliminate some of these problems, such as the choice of where to position the bins and the roughness due to the rectangular bins.  A kernel density estimate based on n observations (Subscript[X, 1]), (Subscript[X, 2]),  ... , (Subscript[X, n]) is defined as

(Subscript[f, h])(
              x) = (1/n) ∑(1/h) K((((Subscript[X, i]) - x)/h))

where the K is a so-called smoothing kernel, typically a smooth probability density such as the standard normal.  The parameter h > 0 controls the level of smoothness:  small values of h produce a rough density estimate, whereas large values produce a very smooth estimate.  The right choice depends on the true underlying density function which is to be estimated, and Silverman (1986) reviews methods that have been proposed for choosing h from data.

My function KernelDensity computes a kernel density estimate, and returns the estimate (Subscript[f, h]) as  a function.  The default options for  KernelDensity allow one to choose the kernel function K (via KernelFunction), the number of points at which the kernel density is evaluated (KernelGrid), and the smoothing parameter h (KernelScale).

Options[KernelDensity]

By default, the kernel scale h is determined automatically by my function KernelScale that assumes the unknown density of the data f is normal with standard deviation  σ and evaluates equation (3.21) from Silverman (1986),

\(TraditionalForm\`h\_opt = \(\(k\_2\^\(\(-2\)/5\)\)(\(∫\(\(K(t)\)\^2\) dt\)\/\(∫f' \(\((t)\)\^2\) dt\))\)\^2\)

where (Subscript[k, 2]) = ∫(t^2) K(t) dtis the second moment of the kernel function.  The resulting smoothing parameter (Subscript[h, opt]) minimizes the mean integrated squared error of the density estimate under the assumption of normality.

For example, a particularly simple kernel function is the biweight.  This function is defined on [-1,1] as K(t) = (15/16) (((1 - (t^2)))^2)and zero elsewhere. (This is the default kernel.)  For this kernel, the optimal smoothing parameter is about  2.78 σ/(n^(1/5)).  Here's the definition.

k[t_] := (15/16) (((1 - (t^2)))^2);
KernelScale[k, -1, 1]

For a Gaussian kernel, the optimal parameter (Subscript[h, opt])is approximately 1.06  σ/(n^(1/5)).

k[
          t_] := ((E^(-((t^2)/2)))/(Sqrt[2 π]));
KernelScale[k, -∞, ∞]

To avoid repeated calculation of the optimal smoothing constant via the ratio of integrals of the kernel, both of these kernels with their optimal parameters are pre-defined.  Further illustration of these calculations in Mathematica appears in Stine (1995).

Returning to the data analysis, KernelDensity by default uses a biweight kernel with (Subscript[h, opt])  defined as above.  The result is an interpolating function on the shown range.  For this example, I've removed the October, 1987 outlier.  I'll use Position to do this to have the indices around for selecting from other series; you can also do this with Select or Cases.

useDates  =  
        Flatten[Position[RETURN,  _?(Abs[#] < 0.05 &)]];

filteredReturn  =  
        RETURN[[useDates]];

kd  =  KernelDensity[filteredReturn]

The symbol kd represents a function which interpolates the kernel density estimate along a default grid of 25 points which are equally spaced on the x-axis.  The domain of the interpolating function is, in this example, the interval [-0.091, 0.092].  This function can be evaluated like any other within its domain of definition.

{kd[0],  kd[ .01],  kd[ .02],  kd[ .03]}

Plots are more interesting.

kplot = Plot[kd[x], {x, -0.06, 0.06}];

The use of a kernel density estimate with the optimal smoothing parameter does not completely avoid the question of how to choose a bin width.  The calculations that lead to the optimal smoothing constant begin with the assumption that the population is normal. Some will feel that the kernel scale chosen by the optimal scheme that defines (Subscript[h, opt]) produces an estimate that is  too ``rough'' in this example. In such situations, one can control the smoothness of the estimate by manually setting KernelScale. Setting the value of this option to 1 and 5 yield rougher and smoother density estimates, respectively.

kdHalf = 
        KernelDensity[filteredReturn, KernelScale ->  .5];

kdTen = 
        KernelDensity[filteredReturn, KernelScale -> 10];

The resulting estimates are rather different.  In the figure that follows, the kernel using (Subscript[h, opt]) is shown as the connected line, with the alternatives using h = 1 in dashed black (its the wiggly one) and h = 5as dashed gray (the normal-looking curve).

kpThree = 
        Plot[{kd[x], kdHalf[x], kdTen[x]}, {x, - .06,  .06}, 
          PlotStyle -> {GrayLevel[0], 
{GrayLevel[0], 
                Dashing[{ .02,  .02}]}, {GrayLevel[ .5], 
                Dashing[{ .02,  .02}]}}
];

Silverman (1986) suggests several more rigorous methods for choosing the scale, but all require further assumptions about the unknown density or substantially more computing.

When overlaid with a histogram, we see that the kernel estimator smooths out irregularities due to the rectangular bins.  Most of its advantage over the histogram as an estimator of the unknown  density comes from removing the rough edges of the histogram's bins.  Both of the estimates, however, indicate a curious ``bump'' on the right side near a return of 0.025.  (The plot "hist" is defined in the previous section on histograms.)

Show[hist, kPlot];

Had we not scaled the histogram to integrate to one, we could not have superimposed these two plots.

Boxplots

We can further augment the plot of the kernel density estimate and histogram by adding a boxplot  to locate the quartiles and outliers.  The latest version of Mathematica now includes a boxplot called BoxWhiskerPlot.  I'll continue to use my own version Boxplot, though.  I prefer the way that I can add it to another plot that I don't find so easy with BoxWhiskerPlot. By default, a boxplot is  horizontal.  (filteredReturn is defined in the previous section.)

Boxplot[filteredReturn];

Because a boxplot summarizes one variable, we need to define how to position it if we want to show it in a graph with two axes.  For this purpose,  BoxCenter is the x-axis or y-axis coordinate that locates the image of the boxplot, and BoxWidth is the width of the box on this same axis.  Since the vertical scale of the histogram of the returns runs from 0 to 35, we locate the box below the x-axis so that it will not obscure the histogram or axis labels and set its width to 5.  It looks the same here, but  now we're set to add the boxplot to the prior histogram and kernel density.

bp = 
        Boxplot[filteredReturn,  BoxCenter ->  -10,  
          BoxWidth ->  5];

The final summary of this subset of the daily returns combines the histogram, kernel density, and boxplot in one display.

Show[hist, kPlot, bp];

Boxplots are particularly useful for comparison.  For example, a sequence of comparison boxplots makes the periods of high variation in the returns clear.   Comparison boxplots are constructed by using a pure function applied to the data binned into 20 cells.  Setting DisplayFunction->Identity avoids having to see the 20 individual boxplots as they are generated.

boxes = (Boxplot[#1[[2]], 
                BoxCenter -> #1[[1]], 
BoxWidth -> 25, BoxVertical -> True, AspectRatio -> 2, 
                DisplayFunction -> Identity] &) /@ 
          Bin[DATE, RETURN, BinPercentage ->  .05];

Show[boxes, PlotRange -> Automatic, 
        DisplayFunction -> $DisplayFunction];

The built-in BoxWhiskerPlot can also show comparison boxplots.

Quantile Plots

As I mentioned before, a weakness of the algorithm that chooses the default smoothing in KernelDensity is that this method assumes that the underlying population is normal.  To see if this assumption is reasonable, we can use a normal quantile plot (or Q-Q plot) of the filtered returns. Mathematica now  includes a version of a quantile.  I need to modify it in order to check for normality, however.

The version of these plots included with Mathematica assumes one is checking to see if the two samples arise from a common population (basically, like a Kolmogorov-Smirnov test).  To get a normal quantile plot, we need the expected normal order statistics.  Those are easy to approximate, and so I've included a modified version of QuantilePlot that allows you to put any of the standard probability distributions as the reference distribution.

The quantile plot includes  a reference line; when the data are a sample from a normal population, the points lie within sampling error of this line.  You can modify this quantile plot to put in quantiles from other standard distributions.  In this exaample, you can see  the data drift away from normality in the tails.

QuantilePlot[filteredReturn];

Although of the upper tail of the data slip away from the reference line at the right of the figure (a region of substantial sampling variation), the distribution of the returns (albeit without October, 1987) appears close to normal.

Scatterplots and Scatterplot Smoothing

Scatterplot Matrix

Although there appears to be little correlation over time in the sequence of returns, one might still  suspect that some correlation exists in the volatility.  The binning function again comes in handy for computing means and standard deviations from the 50 periods of 10 days each.  (This command does not generate  51 bins since it excludes the last incomplete bin of 6 values; the bin size is the floor of the length of the series times the bin percentage.)  The variable useDate is defined in the previous section on kernel density estimates. (If you use the actual returns, the outlier for October dominates the figures.)

bins = 
        Bin[DATE[[useDates]], 
          RETURN[[useDates]], 
          BinPercentage ->  .02];

time = First /@ bins;
m = (Mean[#1[[2]]] &) /@ 
          bins;
s = (StandardDeviation[#1[[2]]] &) /@ 
          bins;

To view these two sequences of summary measures, we would like to examine each versus time as well as plotted versus each other.  My function ScatterplotMatrix generates all three plots, and collects them together in a graphics array with their associated correlations.

ScatterPlotMatrix[m, s, time];

Mathematica now includes its own version of this plot.  Its called PairwiseScatterPlot and it expects the data to be a list of observations, rather than a list of variables.  It gives the redundant form, oriented about the transverse diagonal.

PairwiseScatterPlot[Transpose[{m, s, time}]];

The coincident timing of the outlying values of the mean and standard deviation of the returns is apparent in the last column of two plots which show both series plotted on time.  The plot of the mean versus standard deviation of the returns suggests just how volatile the market became in late 1987.

Scatterplot Smoothing

This section focusses upon functions that describe dependence.  For these examples, the dependence is that in the sequence of 50 standard deviations computed from 10 consecutive trading days.  One period dominates the sequence of standard deviations in the prior scatterplot matrix.  Rather than continue with this outlier, set the extreme value to the average of its neighbors.  The column s is defined in the previous section.

Max[s]

Position[s, Max[s]]

s[[21]] = 
      Mean[s[[{20, 22}]]]

Here is the sequence plot of the ``cleaned'' standard deviations.  With the outlier no longer compressing the rest of the data, this plot shows considerable tracking.  (Compare this plot to that in the lower right corner of the scatterplot matrix.)  ScatterPlot is basically my version of ListPlot, designed to work with data as columns rather than as a list of observations.

ScatterPlot[time, s];

A lag plot gives a better notion of the strength of the correlation.  The scatterplot of $SD_t$ on its lag $SD_{t-1}$ shows a moderate level of period-to-period dependence, confirmed by the correlation estimate which is slightly above 0.5.

lagSD = Drop[s, -1];
SD = Drop[s, 1];

sp = ScatterPlot[lagSD, SD];

Correlation[lagSD, SD]

Whereas correlation measures the strength of linear association as a single value, a scatterplot smoother offers an image of the dependence which does not presume linearity.  My function Smooth, like KernelDensity, builds a function.  The function constructed by Smooth estimates the dependence of the standard deviation at time t on the SD at time t-1. in the form of the conditional expectation E(
          (Subscript[S, t ]) | (Subscript[S, t - 1])). Smooth works by building a sequence of local linear approximations that are joined with a cubic spline.  Near the edges of the data, the smooth is simply the local linear fit.

smth = Smooth[lagSD, SD]

The domain of the resulting function is limited to the  interval indicated in the result that defines the interpolating function.

Plot[smth[x], {x, 0.0064, 0.025},  
        PlotRange -> All];

By default, Smooth constructs each local linear approximation using 40% of the data, as defined by SmoothBinPercentage and is evaluated at about 20 points along the x-axis (the reciprocal of SmoothSkipPercentage).

Options[Smooth]

The wiggles in the smoother, like those in a kernel density which is too rough, are probably random fluctuations.  As with kernel density estimates, the choice is somewhat subjective.  Although computationally intensive methods exist for determining the optimal smoothing parameter, one needs to carefully examine the estimate graphically.  Setting the binning percentage wider to span 60% of the data offers a smoother approximation to the conditional expectation of the local time-dependent standard deviation.

smth60 = Smooth[lagSD, SD, SmoothBinPercentage ->  .6]

I'll plot the second smooth estimate using dashes so that we can superimpose it with the other initial estimate.  With either smoother, the relationship appears a bit nonlinear with the slope tapering off as we move from left to right.

smthPlot = 
        Plot[{smth[x], smth60[x]}, {x, 0.008, 0.025}, 
          PlotStyle -> {GrayLevel[0], 
              Dashing[{0.025, 0.025}]}];

With the data superimposed, the nonlinearity appears small relative to the variation in the sequence of standard deviations.

Show[sp, smthPlot];

Regression

Least Squares

Regression offers a parametric alternative to smoothing.  Mathematica has its own version,  called Regress.  It differs from the more basic program  OLS described here by providing more output and taking longer to run.  One uses OLS in the manner of ScatterPlot, with the covariates or predictors as the first argument, then the response. Two options control the behaviour of OLS.

Options[OLS]

By default,  the result of OLS is a list with the fitted values and residuals.  If the option Residuals is set to  False, the fitted values and residuals are not returned; instead, OLS returns a list of the fitted coefficients.  A contant term is estimated in the linear equation unless the option Constant is set to False.

In this example, the fitted slope of a simple linear regression of $SD_t$ on its lag $SD_{t-1}$ is quite significant.

fitRes = OLS[lagSD, SD];

ScatterPlot accepts the output of OLS as its input and produces the familiar plot of residuals on fitted values.  We can also embellish the residual plot with a smooth function by using the ScatterPlotFunction option and a smooth of the residuals on fitted values.

sRes = Smooth @@ fitRes

ScatterPlot[fitRes, ScatterPlotFunction -> sRes];

The smooth in this residual plot suggests slight deviations from linearity, and once again one might suspect that the slope differs at the extremes of the data. For the curious, we can go further and fit a quadratic polynomial model as a means to capturing the slight curvature. The fit, as one might expect, does not improve significantly (the t-statistic for the added coefficient is small).

OLS[{lagSD, (lagSD^2)}, SD];

We can show this quadratic function with the data as well, once again using the ScatterplotFunction option to superimpose a function on the scatterplot.  Though not significant, the curvature is consistent with the visual impression created by the original nonlinear smoother.

quad[x_] :=  .0047 + 0.84 x - 10.02 (x^2);
ScatterPlot[lagSD, SD, 
        ScatterPlotFunction -> quad];

Robust Regression

While very few methods can accommodate the huge outlier introduced into the analysis by the activity around October, 1987, it is often useful nonetheless to check results of a least squares analysis by using a robust regression.  Street, Carroll, and Ruppert (1988) give a short, readible introduction to the calculations and ideas of robust regression.  A robust regression uses an iterative calculation that performs a sequence of weighted least squares regressions, here starting from an unweighted OLS fit.  Each iteration attempts to downweight outlying values deemed ``unusual'' by the algorithm.  In this example, a robust regression finds little reason to differ from a least squares analysis.

rr = RobustRegression[lagSD, SD];

The results returned by RobustRegression resemble to those of the Mathematica function Regression and are more extensive than those of OLS.  To distinguish among these items, RobustRegression adopts the style of Regression and returns a list of replacement rules.  The labels for each of these replacement rules are shown next.

First /@ rr

To plot the residuals, we have to use rules to extract values as in the following.

robFitRes = {PredictedResponse, 
            FitResiduals} /. rr;

The robust fit is  similar to that from least squares even though several observations have been downweighted by the iteratively reweighted least squares algorithm used to obtain the robust estimates. The plot on the left shows the residuals and fit from the robust analysis, and that on the right is the plot of the robust residuals on the fitting weights used in the final stage of the iterations.

sp1 = 
        ScatterPlot[robFitRes, DisplayFunction -> Identity];
sp2 = 
        ScatterPlot[robFitRes[[2]], 
          Weights /. rr, 
          DisplayFunction -> Identity];

Show[GraphicsArray[{sp1, sp2}], 
        DisplayFunction -> $DisplayFunction];

The curve evident in the plot of the robust residuals on the weights shows how the points become down-weighted as they move away from the fit (represented by the horizontal line at zero in these plots).  The downweighting by default is accomplished by the biweight function b(
          x) = (((1 - (x^2)))^2),  |x|<1,  and zero otherwise.  (The biweight function used in kernel density estimation differs by a multiplicative constant that makes the kernel integrate to one.)  The robust tuning constant is a further option used to scale the residuals defaults to the value giving 95\% efficiency in a Gaussian problem.

Options[RobustRegression]

Further discussion and examples of robust regression in Mathematica appear in Stine (1995).

Discussion

With the proliferation of many specialized packages for statistics, why should one use Mathematica for data analysis?

Few would dispute the claim that Mathematica produces superb graphics that capture the results of a statistical data analysis. Each graph can be easily modified and customized to illustrate each component of the analysis.  The resulting figures are clearly suitable for inclusion in publications.  These same graphics are, unfortunately, rather static.

Until true interactive graphics become standard in Mathematica, however, it will remain the poor cousin of other dedicated statistical computing environments as a foundation for the actual data analysis.  Modern data analysis, as powerfully illustrated by the R-Code software of Cook and Weisburg (1994), relies upon highly interactive graphics such as plot linking, slicing, and brushing.  These methods are not feasible in current releases of Mathematica for most systems.  The methods are feasible in alternative systems, such as the Lisp-Stat  (Tierney  1990, 1995).  Lisp-Stat provides both the capability to link plots so that high-dimensional features of the data are revealed as well as the interface tools, such as slider controls, that control these features.

Nonetheless, the notebook front-end of Mathematica is itself a highly interactive environment.   Notebooks allow one to develop and carefully document the stages of an analysis within Mathematica itself, without having to resort to cutting and pasting results or thoughts on the analysis into other documents.  One should also keep in mind that choices in data analysis are often guided by theory (such as the choice of the kernel smoothing parameter), and the symbolic manipulation features of Mathematica make such calculations routine.

References

Belsley, D.A. (1993).  Econometrics.m:  A package for doing econometrics in \MM.  In {\it Economic and Financial Modeling with Mathematica}, H. Varian, Ed., 300-343.  Springer, New York.

Cook, R. D. and S. Weisburg (1994).  {\it An Introduction to Regression Graphics.}  Wiley, New York.

Silverman, (1986).  {\it Density Estimation for Statistical and Data Analysis}.  Chapman and Hall, London.

Stine, R. A. (1995).  Data analysis using \MM.  {\it Sociological Methods and Research}, 23, 352-372.

Street, J.O., R. J. Carroll, and D. Ruppert (1988).  A note on computing robust regression estimates using iteratively reweighted least squares.  {\it American Statistician}, 42, 152-154.

Tierney, L. (1990). {\it Lisp-Stat}. Wiley, New York.

--- (1995). Data analysis using Lisp-Stat.   {\it Sociological Methods and Research}, 23, 329-351.