Exploratory Data Analysis with Mathematica

EC 313, Spring 2000  C F Baum

We make use of a number of data-analytic tools provided in the DataTools notebook "Data Analysis using Mathematica", by Robert A. Stine, published in Computational Economics and Finance: Modeling and Analysis with Mathematica, Hal Varian, ed., 1996.

This notebook, although originally written for Mathematica 2.2, contains a number of useful facilities for exploratory data analysis. Economic and financial modeling involves developing an appreciation for the features of the data, going beyond summary statistics or simple plots, which enables the modeler to appropriately deal with any unusual aspects of the data and their distribution. In the example illustrated here, we focus on daily stockmarket data from the 1987-1988 period, which contains a distinct outlier: the Crash of 1987, on October 19, which represents a dramatic shift in all characteristics of the data.

The ReadDataFile command expects to find data files containing,on the first line, a description of the data; on the second line, a list of 'variable names'; and on succeeding lines,a value for each of those variables for one observation (in this case, one time period, or trading day).

ReadDataFile["GM8788.DAT"]

Here are the first three and the last of the daily returns observations:

Short[RETURN]

Here are the mean and standard deviation of daily returns, expressed in per cent:

mu = 100 Average[RETURN]

sig = 100 StandardDeviation[RETURN]

Here is the daily stock price; note the dramatic shift in the price series in the October 1987 crash.

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

crash = Flatten[Position[RETURN, Min[RETURN]]][[1]]

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

We may verify that this observation--the 202d--is indeed the crash date, and use that information to label the scatterplot of percentage returns with calendar dates.

calDate[crash]

sp = 
      Scatterplot[100 RETURN, 
        PlotLabel -> "Daily GM Stock Return (%)", 
        Ticks -> {calDate /@ {1, crash, Length[RETURN]}, Automatic}, 
        AxesOrigin -> {0, -25}, AspectRatio -> 0.5]

If returns were normally distributed, we would expect 95% to fall within +/- 2 standard deviations.

Show[sp, 
      GridLines -> {None, {mu - 2 sig, mu + 2 sig}}]

In this case the standard deviation has been inflated by the outliers surrounding the crash, so that somewhat more than 95% fall within the +/- 2 s.d. interval:

withinQ[
        x_] := (mu - 2  sig  ≤  x  ≤  
          mu + 2 sig)

100 N[(Count[100 RETURN, X_?withinQ])/Length[RETURN]]

Let us now focus on monthly returns, using Stine's Bin[] function to group observations which belong to the same month (which we label as 1-24).

byMonth = 
        Bin[N[(YEAR - 1987)*12 + MONTH], 
          100 RETURN];
Length[byMonth]

The byMonth array generates a list of two elements: the bin identifier (month number) and a list of the elements in that bin (in this case the trading-day returns for the month).

First[byMonth]

The summary array is a list of triples: the monthly average, monthly standard deviation, and month identifier.

summary = {Average[#[[2]]], StandardDeviation[#[[2]]], 
              ToString[#[[1]]]} &  /@  byMonth;

summary[[1]]

This plot illustrates the variation in monthly risk and return for each of the 24 months. Month 10 (October 1987) is the clear outlier, even when the several days associated with the crash are averaged into the monthly figures.

TextListPlot[summary, PlotRange -> All, 
      AxesLabel -> {"Avg Return (%)", "SD of Return (%)"}]

If we compute the correlation of risk and return for the data on the plot above, it comes out negative, counter to the predictions of economic theory and, e.g,, the Capital Asset Pricing Model (CAPM): that investors require a higher return to assume greater risk. This would suggest that risk and return should be positively correlated.

avg = Transpose[summary][[1]];
sd = Transpose[summary][[2]];
Correlation[avg, sd]

If we remove the 10th month and consider the correlation of the other 23 months' risk and return values, the correlation has the expected sign.

avg[[10]] *= NA;

Correlation[avg, sd]

We can plot the time series of monthly returns by creating an (approximate)  date variable and using Bin[] to average daily observations over the month:

DATE = DAY + MONTH*30 + (YEAR - 1987)*360;

ListPlot[{#[[1]], Average[#[[2]]]} &  /@  
        Bin[DATE, 100 RETURN, BinWidth -> 30], PlotRange -> All, 
      PlotJoined -> True, 
      PlotLabel -> "Average Monthly Returns on GM (%) for \
1987-1988"]

If we generate a histogram of the daily returns, its scale is dominated by the huge outlier from the Crash:

Histogram[100 RETURN]

Create a filtered return series by considering daily returns less than 10% in absolute value (in this case, only two days' values). This gives a better picture of the distribution of returns.

filtReturn  =  
        Select[100  RETURN, (Abs[#] < 
                10) &];
Length[RETURN] - Length[filtReturn]

hist = 
      Histogram[filtReturn, PlotLabel -> "Filtered GM Return (%)", 
        BinWidth -> 0.25]

We might want to fit a curve through these estimates--an empirical distribution function. It is clear that the filtered returns are not symmetrically distributed, so we would not want to fit a normal curve that imposed symmetry. A kernel density estimate generates a smooth function by averaging 'nearby' observations. Different kernel density functions allow for varying definitions of 'nearby',and the kernel scale controls the smoothness of the resulting function. We compute three kernel density estimates of the filtered returns series: for the default value of the kernel scale parameter, and for two arbitrary values. The higher value of kernel scale generates a smoother function. In all three of the estimated density functions, the 'bulge' at about +2.5% is evident.

kd = KernelDensity[filtReturn]
kd1 = 
      KernelDensity[filtReturn, 
        KernelScale -> 1]
kd5 = KernelDensity[filtReturn, KernelScale -> 5]

kplot = Plot[kd[x], {x, -6, 6}]
kplot15 = 
      Plot[{kd1[x], kd5[x]}, {x, -6, 6}, 
        PlotStyle -> {GrayLevel[0], Dashing[{0.02, 0.02}]}]

A box plot summarizes the distribution of a series, with the box drawn at the 25th and 75th percentile of the series, with the line at the 50th percentile, or median. The horizontal line in this plot has a range of +- 2 standard deviations, with outliers plotted as dots outside that range. Note that these are the filtered returns,excluding the extremes of the Crash.

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

We can combine the histogram and the box plot of the series to illustrate a number of features of the data in a single graph:

Show[hist, bp]

We can also generate a sequence of boxplots, one for each time period, which show in one graph how the median returns and volatility of returns varied over the period.

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

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

Note that these are being generated from the original returns series--not the filtered returns. If we allow the scale to encompass the unusual outcomes of October 1987, the other months' variations become trivial:

Another summary measure of a univariate distribution is the Quantile-Quantile (QQ) Plot. With the quantiles of the Normal distribution as a benchmark, we may consider how the distribution of filtered returns compares with the Normal.Small variations are observed in the extreme tails.

QQPlot[filtReturn]

Contrast this with the distribution of the unfiltered returns--which only contains two additional observations! The returns series seems to diverge from the Normal to a much greater degree.

QQPlot[RETURN]

Let us consider how the standard deviation of returns, or volatility, varied over the sample. First, with respect to the original returns series:

bins = 
        Bin[DATE, RETURN, 
          BinPercentage -> 0.02];
time = First  /@ bins;
avg =  Average[#[[2]]] &  /@  
          bins;
stddev = StandardDeviation[#[[2]]] &  /@  
          bins;

Scatterplot[time, stddev, PlotJoined -> True, 
      AxesOrigin -> {0, 0}]

We may also display the three series' comovements against each other (as well as their bivariate correlations) in a scatterplot matrix:

ScatterplotMatrix[avg,stddev,time]

If we now locate the volatility outlier, corresponding to the Crash,and average its value with the adjacent observations, we get quite a different picture of the volatility of returns:

Position[stddev, Max[stddev]][[1]]

stddev[[21]] = Average[stddev[[{20, 22}]]]

Scatterplot[time, stddev, PlotJoined -> True, 
      AxesOrigin -> {0, 0.007}]

Let us now consider how volatility may be bunched, or clustered, with high  (or low) volatility persisting. Using this filtered volatility series, compute a lagged value, and plot the current value vs. the lagged value:

lagSD = Drop[stddev, -1];
SD = Drop[stddev, 1];
sp = Scatterplot[lagSD, SD, AspectRatio -> 1]

We may calculate the correlation between successive periods' volatility, and find it to be substantial:

Correlation[SD, lagSD]

Correlation examines the linear association of two series. A more general concept is that produced by a scatterplot smoother, which like a kernel density estimate generates a function relating the two series from a number of 'nearby' values, without the constraint of linearity.  If we plot this smoothed function (superimposed on the scatterplot from which it was generated), we can detect a mild nonlinearity in the relationship between contemporaneous and lagged volatility:

smth = Smooth[lagSD, SD]

smthPlot = 
        Plot[smth[x], {x, 0.008, 0.025}, AspectRatio -> 1, 
          DisplayFunction -> Identity];

Show[sp, smthPlot]

We may also consider this relationship between current and lagged volatility in a regression context, and generate a smoothed scatterplot of the regression residuals (or in-sample prediction errors). Again, a mild nonlinearity is apparent, since the smoothed line would be horizontal at zero if the distribution of the errors conformed to that of white noise.

fitRes = OLS[lagSD, SD];

sRes = Smooth  @@  fitRes

Scatterplot[fitRes, ScatterplotFunction -> sRes]

Robust regression, also included in Stine's package, could be used to estimate a relationship for which the data contained extreme values, such as those of the Crash of 1987.