Morgan Stanley
Quantitative and Mathematical Finance Conference
Prague October 17th - 18th, 2005

Selected Statistical Methods with Applications to Finance

Jan Hurt *)

Charles University of Prague, Czech Republic

*) The research has been supported by Grant of Czech Ministry of Education MSM 0021620839 and by Morgan Stanley.

1 Introduction

This contribution is devoted to some of the multivariate statistical methods and comments to its use in the field of finance. Some of them are already routinely used in finance, the application of the others is not apparent. Nevertheless, it is a strong  author's belief that all of the methods have a great potential in explanation and prediction of many financial features. By the author's opinion, the discriminant analysis is nowadays the most powerful and useful tool among all multivariate methods mentioned here. An interesting method with possible use is multidimensional scaling which will be mentioned as the last one. 

Beyond a pile of standard software devoted to the methods of multivariate statistical analysis there are also systems which are called Computer Algebra Systems (CAS). With help of a qualified statistician, these systems may highly contribute to the deep analysis of financial  phenomena due to their very simple and powerful computational environment and interactive approach. We will illustrate some of the methods below via CAS Mathematica, the system of doing mathematics by computer created by Stephen Wolfram.

1.1 General preliminaries

The leitmotive of all multivariate statistical methods is the reduction of information contained in the data. Here we suppose that the data are represented by n observations (Subscript[X, 1]),  ... , (Subscript[X, n]) of a p-dimensional vector X. Most of the methods are based on studying the covariance structure of the data.  

Vector X consists of measurable characteristics of individuals (firms, clients) in question. A typical example is credit risk modelling and analysis. Here the observable variables (the components of vector X) may be turnovers, equity, size of credit, rating, branch of enterprising, classified credit, recovery rate, etc.

Despite almost all of the methods are dealing with known characteristics of random phenomena in question, in practise the user is forced to use some appropriate estimates of unknown characteristics, in case of known analytical form of the distributions applies to the unknown parameters.

1.2  Mathematical preliminaries

It may be considered as a very naive approach but a crucial role of the analysis plays the so called spectral decomposition of a symmetric real (p x p) matrix A. In such a case there exist an orthogonal matrix P and a diagonal matrix Λ such that

A = P Λ (P^T).

On the diagonal of Λ we have the eigenvalues (Subscript[λ, 1])≥ ...≥(Subscript[λ, p]) of A and the columns of P are the corresponding orthonormal eigenvectors (Subscript[x, 1]),  ... , (Subscript[x, p]). An alternative form is then

A = (∑_(i = 1)^p) (Subscript[λ, i]) (Subscript[x, i]) 
            ((Subscript[x, i])^T).

Note, that if the eigenvalues are not all different, the previous decomposition is not unique. All the methods described below do not depend on the possible nonuniqueness, however.

2 Principal Components

The objective of Principal Components Analysis (PCA) is to express the statistical information contained in random vector X by k ≤ p linear combinations of this random vector with at least as possible loss of information. The question is how large should k be and how to construct the mentioned linear combinations to recover the random vector X. 

Suppose that the random vector in question possesses covariance matrix Σ (necessarily positive semidefinite) with eigenvalues (Subscript[λ, 1])≥ \(TraditionalForm\`\(\(...\)\(≥\)\(λ\_p\)\(≥\)\(0\)\)\), and the  corresponding orthonormal eigenvectors (Subscript[x, 1]),  ... , (Subscript[x, p]) so that the covariance matrix Σ may be expressed as

\(TraditionalForm\`Σ = ∑\+\(i = 1\)\%p\( λ\_i\) \(x\_i\) x\_i\^T\).

Random variable (Subscript[Y, i]) = ((Subscript[x, i])^T) X is called i-th principal component, i = 1,  ... , p. The variance of the i-th principal component is var (Subscript[Y, i]) = (Subscript[λ, i]). If the rank of the covariance matrix is r then the last p - r principal components are constants almost surely. The total variance of random vector X is

(σ^2) = (∑_(i = 1)^p) 
              var (Subscript[X, i]) = (∑_(i = 1)^p) 
                var (Subscript[Y, i]) = (∑_(i = 1)^p) (Subscript[λ, i]).

Surely, (σ^2) is a certain measure of the variability of X. The decomposition shows how large portion of the overall variability is explained by the i-th principal component. Usual practical approach consists in construction of the first s components, say, and looking how large part of the total variance is expressed by principal components (Subscript[Y, 1]),  ... , (Subscript[Y, s]) by making use of the ratio (∑_(i = 1)^s) (Subscript[λ, i])/(σ^2). If  the last expression is sufficiently close to 1, the first s principal components are considered to be sufficient to explain to total variability. It is surprising that in properly formed models, only a few components (three, say), explain over 90 per cent of the total variability.

3 Discriminant Analysis

3.1 Introduction

Despite discriminant analysis (DA) is treated almost in any textbook a little emphasis is put on the fact that DA in principle, is a two-stage procedure based on multivariate observations of a (p x 1) random vector X. Even in some specialized program packages we can sometimes see only the first stage. But the second stage is the most important for decision making!

At the first stage we have so called training set of n individuals, say, represented by (p x 1) observations (Subscript[X, 1]), ... (Subscript[X, n]). In what follows let us identify the individuals 1,  ... , n with their respective observations. At this stage, we have also the information to which of k populations (classes, groups) an individual belongs to. The populations will be marked 1,  ... , k. The information about belonging of an individual (Subscript[X, j]) to a population is expressed by the indicator (tag) (Subscript[i, j]), (Subscript[i, j])  ∈  {1,  ... , k}, j = 1,  ... , n. The input data for the first stage are in the form data={{(Subscript[X, 1]), (Subscript[i, 1])}, ... {(Subscript[X, n]), (Subscript[i, n])}}, a three-dimensional structure. The purpose of the first stage is to build up criteria for classifying other individuals into respective groups in the second stage. The methods considered here are either purely nonparametric (Mahalanobis), semiparametric (Bayesian), or parametric ones based on the assumption of normality of the data involved. In the latter case, linear or quadratic DA is under consideration. 

In practice the exact distributions are rarely known. In better case we know just an analytical form of the distributions involved. In case we know the analytical form of the distributions (parametric case, we have to estimate the unknown parameters. In this case we deal with the multinormal distribution. The populations may differ either in means or both in means and  covariance structure. In other cases we use the Bayesian or Mahalanobis approach.

At the second stage, classification, we use discriminant scores created at the first stage to classify newly incoming individuals with known observations X (and apparently unknown class) into one of the k classes. Usually an individual is assigned to population with maximum or minimum discriminant score.

Finally, it is of interest to evaluate the performance of the discrimination rules. The simplest method is the resubstitution method. We allocate the individuals from the training set using the discrimination rules and look how often the data are misclassified. The drawback of this method is that the individual is classified by rules which have been created utilizing the information of that individual. More sophisticated methods are know as cross-validation, jackknife, (the last mentioned we consider as a very satisfactory approach) see [1], [2], e. g.

In what follows let (Subscript[L, i])(x) denote the likelihood or loglikelihood function of the population i, i = 1,  ... , k. These functions help us to construct so called discriminant scores which are used in the second stage. The common form of discriminant scores is (Subscript[S, 1])(x),  ... , (Subscript[S, k])(x) and the newly incoming individual x is assigned to population arg (max_i) (Subscript[S, i])(x).

3.2 Preliminaries

3.2.1 Test data

In Mathematica environment we create three two-dimensional populations of sizes 40, 50, 60, respectively. We can consider the three populations as the three industries. The two-dimensional observations may be the maximum and minimum expected rate of return during a given period. The respective risk is given by the covariance matrices.

Needs["Statistics`"];
Needs["\
Statistics`MultiDescriptiveStatistics`"];
SeedRandom[13];
μ1 = {5.006, 
        3.428}; μ2 = {5.936, 2.770}; μ3 = {6.588, 
        2.974};
Σ1=(({{0.1218, 0.0972}, {0.0972, 0.1408}}));
Σ2=(({{0.2611, 0.0835}, {0.0835, 0.0965}}));Σ3=(({{0.3963, 0.0919}, {0.0919, 0.1019}}));
d1 = 
        MultinormalDistribution[μ1, Σ1];
d2 = 
        MultinormalDistribution[μ2, Σ2];
d3 = 
        MultinormalDistribution[μ3, Σ3];
x1 = 
        RandomArray[d1, {40}];
x2 = 
      RandomArray[d2, {50}]; x3 = RandomArray[d3, {60}];

3.2.2 A graphical representation using Help function of  FindClusters

Note that starting from Mathematica version 5.1 the user can use an excellent Add-On function FindClusters. We will mention it later in Section on Cluster analysis. Here we only use the graphical environment of it.

ClusterPlot[cl_, opts___] := 
        Module[{len, colors = {Red, Blue, Green, Purple, Cyan}}, 
          colors = 
            PadRight[colors, Length[cl], colors]; 
Show[
            Graphics[{AbsolutePointSize[5], 
                Transpose[{colors, Map[Thread[Point /@ #] &, cl]}]}, opts, 
              Frame -> True]]];

ClusterPlot[{x1, x2, x3}];

3.2.3 Data equipped with tags

The indicators of the populations are added:

newdata = 
        Join[Map[Join[#, {1}] &, x1], Map[Join[#, {2}] &, x2], 
          Map[Join[#, {3}] &, x3]];

newdata

3.3 First stage

3.3.1 Calculations

Description of the symbols

n ... the size of the training set
populationmarks ... tags of populations
k ... number of populations
populations ... the data of the training set partioned into k classes together with their tags
populationsred ... the same structure as above but without tags (reduced data); necessary for estimating unknown parameters
frequencies ... the sizes (Subscript[n, 1]), ... (Subscript[n, k]) of populations in the training set
priorrel ... prior probabilities of the populations bases on relative frequencies in the training set
covmatrices ... covariance matrices of the populations (unbiased estimates)
pooledcovariance matrix ... the unbiased estimate of the covariance matrix of the whole data assuming the equal covariance matrices in the populations

\(\(\(Needs["Statistics`MultiDescriptiveStatistics`"]\)\(
\)\( (*\ \(data = {{11. , 12, 1}, {33, 332, 3}, {22, 222, 
              2}, {333, 3332, 3}, {3333, 33332, 3}, {111, 1112, 1}, {222, 
              2222, 2}, {1111, 11112, 1}, {33333, 333332, 
              3}};\)*) \)\(
\)
    \) (*\ Size\ of\ the\ training\ set\ *) \)
n = Length[newdata]
populationmarks = 
      Union[Last[Transpose[newdata]]]

      
k = Length[populationmarks]

      
populations = 
        Map[Cases[newdata, {__, #}] &, 
          populationmarks];
  
frequencies = Map[Length, populations]

      
priorrel = frequencies/n

      
populationsred = 
        Drop[populations, {}, {}, {-1}] // 
          N;
means = Map[Mean, populationsred];
covmatrices = 
        Map[CovarianceMatrix, populationsred];
invcovs = Map[Inverse, covmatrices];

      
pooledcovariancematrix = 
        Apply[Plus, (frequencies - 1)*covmatrices]/(n - 
              k);

means

Map[Mean, populationsred]

covmatrices

Map[CovarianceMatrix, populationsred]

3.3.2 Procedures

Linear DA

Clear[lDA];
lDA[populationsred_, pooledcovariancematrix_, prior_]:=
	\(Module[{means, invcov, 
              scores}, 
		means = 
              Map[Mean, populationsred]; 
		invcov = 
              Inverse[pooledcovariancematrix]; 
		scores = {means . invcov, \(-0.5\)\ Map[# . invcov . # &, means] + 
                  Log[prior]}
		]\)

Quadratic DA

Clear[qDA];
\(qDA[populationsred_, covmatrices_, prior_] := 
          Module[{means, invcovs, dets, scores}, 
means = 
              Map[Mean, populationsred]; 
	
            invcovs = Map[Inverse, covmatrices]; 
	
            dets = Map[Det, covmatrices]; 
	
            scores = {\(-\(1\/2\)\) Log[dets] + Log[prior], means, 
                invcovs}
	\ ]\)

Mahalanobis DA

Clear[mahalDA];
\(mahalDA[populationsred_, 
          pooledcovariancematrix_] := 
        Module[{means, invcov, scores}, 
means = 
            Map[Mean, populationsred]; 
	
          invcov = Map[Inverse, pooledcovariancematrix]; 
	
          scores = {means, 
              invcov}
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ];\)

Maximum Likelihood DA

Clear[mlDA];
\(mlDA[populationsred_, 
          covmatrices_] := 
Module[{means, invcovs, dets, 
            scores}, 
	
          means = Map[Mean, populationsred]; 
	
          invcovs = Map[Inverse, covmatrices]; 
	
          dets = Map[Det, covmatrices]; 
	
          scores = {dets, means, 
              invcovs}
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ];\)

Bayesian DA

Clear[bayesDA];
\(bayesDA[populationsred_, covmatrices_, 
              prior_] := 
            Module[{means, invcovs, dets, scores}, 
	
              means = Map[Mean, populationsred]; 
	
              invcovs = Map[Inverse, covmatrices]; 
	
              dets = Map[Det, covmatrices]; 
	
              scores = {\(-\(1\/2\)\) Log[dets] + Log[prior], means, 
                  invcovs}\ ];\)

lDA[populationsred, pooledcovariancematrix, priorrel]

3.3.3 Construction of the criteria

Possible options:
priors->uniform, relativeFrequencies (default)
method-> linearDA (default),quadraticDA, MahalanobisDA, MLDA, BayesianDA

Clear[criteriaDA];
Options[criteriaDA] = {method -> linearDA, 
            priors -> relativeFrequencies};
  
\(criteriaDA[data_, 
        opts___] := 
Module[{n, k, populationmarks, 
          populations, frequencies, priorrel, covmatrices, 
          pooledcovariancematrix, met, pri, 
          prior (*\(,\)\(lDA\)\(,\)\(qDA\)\(,\)\(mahalDA\)\(,\)\(mlDA\)\(,\)\(bayesDA\)*) }, 
n = Length[data]; \  (*\ 
          the\ size\ of\ the\ training\ set\ *) 
populationmarks = 
          Union[Last[Transpose[data]]]; 
k = 
          Length[populationmarks];  (*\ 
          No . \ of\ populations\ *) 
populations = 
          Map[Cases[data, {__, #}] &, 
            populationmarks]; 
frequencies = 
          Map[Length, populations]; 
priorrel = 
          frequencies/n; 
 (*\ 
          Removes\ the\ marks\ of\ populations\ from\ the\ data\ set\ but\ saving\ the\ structure\ of\ populations\ *) 
populationsred = Drop[populations, {}, {}, {\(-1\)}] // N; 
covmatrices = 
          Map[CovarianceMatrix, populationsred]; 
 (*\ 
          Creates\ pooled\ covariance\ matrix\ *) 
pooledcovariancematrix = 
          Apply[Plus, \((frequencies - 1)\)*covmatrices]/\((n - k)\); {met, 
            pri} = \({method, priors} /. {opts}\) /. 
            Options[criteriaDA]; 
If[pri === uniform, 
          prior = 1/k, prior = priorrel]; 
Switch[
          met, 
	\ \ linearDA, \ {met, 
              lDA[populationsred, pooledcovariancematrix, prior]} // 
            Return, 
	\ quadraticDA, {met, 
              qDA[populationsred, covmatrices, prior]} // 
            Return, 
	\ MahalanobisDA, {met, 
              mahalDA[populationsred, covmatrices]} // 
            Return, 
	\ MLDA, {met, 
              mlDA[populationsred, covmatrices]} // 
            Return, 
	\ BayesianDA, \ {met, 
              bayesDA[populationsred, covmatrices, prior]} // 
            Return]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\)

crlin = criteriaDA[newdata]

crqua = criteriaDA[newdata, method -> quadraticDA]

crmahal = criteriaDA[newdata, method ->  MahalanobisDA]

crML = criteriaDA[newdata, method -> MLDA]

crbayes=criteriaDA[newdata, 
        method -> BayesianDA]

3.4 Second stage

3.4.1 Construction of classifiers

Important note: The second argument of the following procedures should be criteriaDA[[2]]

Linear DA

Clear[critlDA];
\(critlDA[individual_, 
        crit_] := 
	
      Module[{scores, ma}, 
	
        scores = 
          crit[[1]] . individual + 
            crit[[2]]; 
	
        ma = Max[
            scores]; 
	\((Position[scores, ma] // 
              Flatten)\)[[1]]
		\ \ \ \ ]\)

critlDA[{15, 5}, 
      crlin[[2]]]

Map[critlDA[#, crlin[[2]]] &, 
      x1]

Quadratic DA

Clear[critqDA];
\(critqDA[individual_, 
        crit_] := 
	
      Module[{scores, ma, means, invcovs, au, k}, 
	
        means = crit[[2]]; 
	k = Length[means]; 
	
        invcovs = 
          crit[[3]]; 
	au = Map[individual - # &, means]; 
	
        scores = 
          crit[[1]] - 
            0.5  Map[
                au[[#]] . 
                    invcovs[[#]] . 
                    au[[#]] &, 
                Range[k]]; 
	
        ma = Max[
            scores]; 
	\((Position[scores, ma] // 
              Flatten)\)[[1]]
		\ \ \ \ ]\)

critqDA[{5, 3}, 
      crqua[[2]]]

Map[critqDA[#, crqua[[2]]] &, 
      x1]

Mahalanobis DA

Clear[critmahalDA];
\(critmahalDA[individual_, 
        crit_] := 
	
      Module[{scores, ma, means, invcovs, au, k}, 
	
        means = crit[[2]]; 
	k = Length[means]; 
	
        invcovs = 
          crit[[3]]; 
	au = Map[individual - # &, means]; 
	
        scores = \(-0.5\) 
            Map[au[[#]] . 
                  invcovs[[#]] . 
                  au[[#]] &, 
              Range[k]]; 
	
        ma = Max[
            scores]; 
	\((Position[scores, ma] // 
              Flatten)\)[[1]]
		\ \ \ \ ]\)

Map[critmahalDA[#, crqua[[2]]] &, 
      x2]

Maximum Likelihood DA

Clear[critMLDA];
critMLDA[individual_, crit_] := 
      critqDA[individual, crit]

Map[critMLDA[#, crqua[[2]]] &, 
      x3]

Bayesian DA

Clear[critbayesDA];
\(critbayesDA[individual_, crit_] := 
      Module[{scores, ma, means, invcovs, au, k}, 
	
        means = crit[[2]]; 
	k = Length[means]; 
	
        invcovs = 
          crit[[3]]; 
	au = Map[individual - # &, means]; 
	
        scores = 
          crit[[1]] - 
            0.5  Map[
                au[[#]] . 
                    invcovs[[#]] . 
                    au[[#]] &, 
                Range[k]]; 
	
        ma = Max[
            scores]; 
	\((Position[scores, ma] // 
              Flatten)\)[[1]]
		\ \ \ \ ]\)

Map[critbayesDA[#, crqua[[2]]] &, 
      x3]

3.4.2 Procedure for classification

Variable individual is the vector of observations on an individual. Variable criteria should be obtained using procedure criteriaDA.

Clear[classifyDA];
classifyDA[individual_, criteria_]:=
	\(Module[{met, crit}, 
	
        Switch[criteria[[1]], 
	linearDA, 
          critlDA[individual, 
              criteria[[2]]] // 
            Return, 
	quadraticDA, 
          critqDA[individual, 
              criteria[[2]]] // 
            Return, 
	mahalDA, 
          critmahalDA[individual, 
              criteria[[2]]] // 
            Return, 
	MLDA, 
          critMLDA[individual, 
              criteria[[2]]] // 
            Return, 
	bayesDA, 
          critbayesDA[individual, 
              criteria[[2]]] // 
            Return
	\ \ \ \ \ \ \ \ \ \ \ ]
	\ \ \ \ \ \ \ \ \ \ \ ]\)

Map[classifyDA[#, crlin] &, x1]

3.5 Misclassification

3.5.1 Resubstitution method

allocationresults = 
      Table[Map[classifyDA[#, crlin] &, 
          populationsred[[i]]], {i, 
          Length[populationsred]}]

Probabilities of misclassification

probabilitiesofmisclassification = 
      1 - Table[
          Count[allocationresults[[i]], 
              i]/Length[
              allocationresults[[
                i]]], {i, 
            Length[allocationresults]}]

3.5.2 Cross-Validation

As mentioned earlier, there are problems with correct interpretation of the misclassified items. Obviously, the resubstitution method, by its essence, favours the so called pathological cases since the decision criteria include all of the objects into account and even the moderate  have an extreme influence on decision making. Among the various methods how to overcome the problem, jackknife (in fact, some 59 years old idea seemingly coming from Quenouille) seems to be the most correct one for an empirical testing the efficiency of DA.

4 Factor Analysis

4.1 Introduction

Factor Analysis (FA) is particularly popular in a non-mathematical society (human science, sociology, medicine). There is an authors' opinion that the use of FA requires a good understanding of the subject by the two parties, the researcher and the statistician. The examples in textbooks are promising enough to persuade almost any researcher or financial analyst without thorough practise in statistics about the usefulness of FA and about its capability to recognize hidden dependencies and resolve any problem. In case of FA, the author is a bit pessimistic because looking on the results of calculations, hardly the researcher (a financial analyst in our case) or the statistician are able to give the interpretation of the common factors even after many sophisticated rotations, one of the principal method that should help to better interpretation of the extracted factors. Nevertheless, the brief outlook on the methodology follows. Note that almost all methods of FA are based on the analysis of the covariance structure.

4.2 Factor Model

From the mathematical point of view, FA is a generalization of the model of principal components. As in PCA, so called common factors are considered, but additionally, specific factors for each variable of X are included into the model. The leading idea is that every particular variable of X can be expressed as a linear combination of common factors plus the value of the specific variable, in this scheme called specific scores.  The specific factors make the difference between PCA and FA, and apparently substantionally complicate both mathematical calculations and the interpretation of the results. These ideas will now be translated into the terms of covariance structure. Without any loss of generality let us suppose that the data are properly standardized so that instead of the covariance matrix of the data we will analyze the corresponding correlation matrix, R, say. In finance R is often the correlation matrix of the of returns. The correlation structure may be expressed as in the form of decomposition

\(TraditionalForm\`R = Λ\ Λ\^T + Ψ\).

Here Λ is called the matrix of factor loadings and  Ψ is a diagonal matrix with the variances of specific factors on it. The diagonal elements  ((Subscript[ψ, i])^2) of Ψ are called specificities while the ((Subscript[h, i])^2) = 1 - (ψ^2) are called communalities. 

The objective of FA (from the mathematical point of view) is (i) to estimate the communalities and the loadings, (ii) to give an interpretation to the obtained. The step (i) is now usually done by modified method of principal components. If some pieces of information are known on the distributions involved, then the principle of maximum likelihood is usually adopted.

Since sometimes (the usual phrase in a textbook on FA) or rarely (the author's opinion) the proper interpretation of the factors may be realized, some finer techniques are used. Among tens of them, the most popular are rotations. However, we can hardly assure that a successful result will be attached.

4.3 Illustration to Finance: APT as a Factor Model

4.3.1 Introduction

In various textbooks the reader may find the so called "necessary" assumptions for developing the Markowitz minimum-variance portfolio, Arbitrage Pricing Theory (APT), etc. Most of these assumptions only help the economists to understand the topic but for the development of a mathematical model or even theory are nuisance. See [4] for a discussion about the assumptions on efficient markets, e. g. The mathematical model has nothing to do with the verbal assumptions.s

4.3.2 Factor Model of APT

After what has been said above, it is a simple task to formulate the model of FA in APT environment. The idea is  consider some imaginary common factors that have an influence on the returns. It is therefore on people engaged in this research to do perhaps the most difficult part. To find what is the meaning of the extracted factors.

4.3.3 Regression Model of APT

In this case, the variables and their meaning is known. It is a direct generalization of the Capital Asset Pricing Model (CAPM). Let us recall that in CAPM it is supposed that the rate of return of an asset depends in a linear way on the market rate of return only.  In the generalized form of CAPM, i. e. APT, the rates of returns of various firms are simply considered as (possibly) linear functions of some concrete variables like inflation, growth in industrial production, change in unexpected inflation, risk premium as the difference in yields of corporate bonds and long-term Treasury bonds, and as in CAPM rate of return on market portfolio measured by NYSE index, say.

5 Cluster Analysis

5.1 Introduction

In fact, Cluster Analysis (CA) is by many statisticians hardly considered as a statistical method. Nevertheless, it may help a much in the process of investigation of a mess in multivariate data. CA is something between modern data mining and the classical Sir Ronald A. Fisher's statistical approach. As mentioned elsewhere, Mathematica starting version 5.1 has a pretty nice function FindClusters with a perfect documentation and some very  useful procedures which do not appear even in qualified texts. The detailed information on demand.

5.2 Hierachical Models

5.3 Nonhierarchical Models

6 Multidimensional Scaling

6.1 Introduction

The purpose of the multidimensional scaling (MS) is in some respect reverse to usual statistical methods. Rather than observations we have the distances or, better to say, dissimilarities among objects in a real linear space of possibly unknown dimension. The purpose of MS is to create a configuration of points in a linear space such that the actual distances between points are as close as possible to the prescribed ones. It is shown that if the distance matrix is Euclidean then such a configuration exists and since the proof is in fact the algorithm for constructing the desired configuration. The construction uses the spectral decomposition of a symmetric matrix. In Mathematica this is done by use of functions Eigenvalues and Eigenvectors.

Notation: n number of points, p 

makeconfiguration::usage="makeconfiguration[dit_,  pdim___] gives a configuration of points in a p-dimensional  Euclidean space where p is the number of positive eigenvalues of the respective matrix, or, if pdim given, p:=pdim. The construction uses the distances given in the form of lower diagonal list dit (without zeros on the diagonal), i.e., dit =  {{(Subscript[d, 21])}, {(Subscript[d, 31]), (Subscript[d, 32])}, {(Subscript[d, 41]), (Subscript[d, 42]), (Subscript[d, 43])}, ... }. The configuration is in the form of (nxp) matrix with the centre of gravity (the average) 0."

newdistancematrix::usage="newdistancematrix[pot_] returns the distance matrix of the points given in the (nxp) matrix pot. pot should be of the same form as the output of makeconfiguration."

newdistancediag::usage="newdistancediag[pot] returns the lower diagonal of the distance matrix of the points given in the (nxp) matrix pot. pot should be of the same form as the output of makeconfiguration."

stress::usage="stress[distancediag, distancehatdiag] returns     
ψ = \(TraditionalForm\`∑\+\(i = 1\)\%n\(∑\+\(j = 1\)\%n\((d\_ij\)\)\) - \(TraditionalForm\`d\&^\)\(TraditionalForm\`\_ij\))^2, also called raw stress, where (Subscript[d, ij]) are the original distances and  \(TraditionalForm\`d\&^\)\(TraditionalForm\`\_ij\) are the "estimated" or reproduced distances from the configuration in p dimensions. Both of the parameters should be in the form of lower diagonal structures {{(Subscript[d, 21])}, {(Subscript[d, 31]), (Subscript[d, 32])}, {(Subscript[d, 41]), (Subscript[d, 42]), (Subscript[d, 43])}, ... } (without zeros) of the same dimension.

6.2 Implementation

6.2.1 Illustration of the method -  road distances between Czech cities

Praha (Pr), Podebrady (Po), Hradec Kralove (Hr), Jihlava (Ji)
The distances {{Po-Pr},{Hr-Pr,Hr-Po},{Ji-Pr,Ji-Po,Ji-Hr}}:

dit = {{54}, {112, 58}, {123, 102, 110}};
  
n = Length[dit] + 1;
lowdiag = 
        Join[{Table[0, {n}]}, 
          Table[Join[dit[[i]], 
              Table[0, {n - 
                    Length[dit[[
                        i]]]}]], {i, 
              n - 1}]];

      
Print["The distance matrix"]
(d = lowdiag + Transpose[lowdiag]) // MatrixForm
\( (*\ nxn\ matrix\ of\ ones\ *) \(j = 
        IdentityMatrix[n] - IdentityMatrix[n] + 1;\)
 (*\ 
      centering\ matrix\ *) \)
h = IdentityMatrix[n] - j/n;
a = -d*d/2. ;
b = h . a . h;
Print["Eigenvalues:   ", 
      lambda = Eigenvalues[b]  ]
Print["Number of positive eigenvalues:  ", 
      p = Count[lambda, x_ /; x > 0]]
charvect = Take[Eigenvectors[b], p];
points = 
        Take[Sqrt[lambda], p]*charvect  //  Chop;

      
Print["Points:  ", (pot = Transpose[points]) // 
        MatrixForm]
Print["Actual distances between points  :", Table[
              Apply[Plus, (pot[[r]] - 
                      pot[[
                        s]])^2], {r, n}, {s, n}]  //  
            Sqrt  //  Chop // MatrixForm]

Here is the plot of points:

ListPlot[pot, Prolog -> PointSize[0.04]];

6.2.2 The procedures

Clear[makeconfiguration];
makeconfiguration[dit_, pdim___] := 
Module[{n, lowdiag, d, j, h, a, b, lambda, 
     p, charvect, points, pot, i, r, s}, 
     n = Length[dit] + 1; 
     lowdiag = Join[{Table[0, {n}]}, Table[Join[dit[[i]], 
         Table[0, {n - Length[dit[[i]]]}]], {i, n - 1}]]; 
     d = lowdiag + Transpose[lowdiag]; 
     j = IdentityMatrix[n] - IdentityMatrix[n] + 1; 
     h = IdentityMatrix[n] - j/n; 
     a = (-d)*(d/2.); 
     b = h.a.h; 
     lambda = Chop[Eigenvalues[b]]; 
     p = Count[lambda, x_ /; x > 0]; 
     p = If[pdim === {}, p, pdim]; 
     Print["Dimension of the space:  ",p];
     charvect = Take[Eigenvectors[b], p]; 
     points = Chop[Take[Sqrt[lambda], p]*charvect]//Chop; 
     pot = Transpose[points]]
     
Clear[newdistancematrix];
newdistancematrix[pot_] := 
Module[{r, s, n}, 
    	n = Length[pot]; 
        Chop[Sqrt[Table[Plus @@ ((pot[[r]] - pot[[s]])^2), {r, Length[pot]}, 
        {s, Length[pot]}]]]]
        
Clear[newdistancediag];
newdistancediag[pot_] := 
Module[{r, s, n}, n = Length[pot]; 
   Chop[Sqrt[Table[Plus @@ ((pot[[r]] - pot[[s]])^2), {r, 2, n}, 
        {s, 1, r - 1}]]]]
        
Clear[stress];
stress[distancediag_, distancehatdiag_] := 
   2*Plus @@ Flatten[(distancediag - distancehatdiag)^2]

6.3.3 Example (city distances again)

dit1 = {{54}, {112, 58}, {123, 102, 
            110}};
(pot1 = makeconfiguration[dit1]) // 
      MatrixForm
(dhat1 = newdistancediag[pot1]) // 
      MatrixForm
stress[dit1, dhat1]

6.3.4 Example (200 random points in the five-dimensional space)

randompoints[n_  , k_  ] := 
      Table[Table[Random[], {k}], {n}]

Prescribed dimensions

SeedRandom[13];
dit3 = 
        newdistancediag[randompoints[200, 5]];
\(Table[
	
      Print[Timing[\(pot3 = 
              makeconfiguration[dit3, 
                j];\)]]; 
\ \ \ \ \ \ \ \ dhat3 = 
        newdistancediag[pot3]; 
\ \ \ \ \ \ \ \ stress[
        dit3, dhat3], {j, 2, 5}]\)

Free choice of the dimension

SeedRandom[13];
dit3 = 
        newdistancediag[randompoints[200, 5]];
Timing[pot3 = makeconfiguration[dit3];]
dhat3 = newdistancediag[pot3];
stress[dit3, dhat3]

7 Conclusions

The DA procedures are restricted either to purely nonparametric or normal case. In maximum likelihood and Bayesian DA, a natural generalization of procedures is in introducing the other multivariate distributions as parameters. A possible candidate for that is a class of copulas.

8 References

[1] Christensen, R.: Linear Models for Multivariate, Time Series, and Spatial Data. Springer-Verlag. New York, 1991.
[2] Mardia, K. V., Kent, J. T., Bibby, J. M.: Multivariate Analysis. Academic Press. London, 1979.
[3] Rao, C. R.: Linear Statistical Inference and Its Applications. 2nd edition. Wiley. New York, 1973.
[4] Dupacova, J., Hurt, J. Stepan, J.: Stochastic Modeling in Economics and Finance. Kluwer Academic Publishers. Dordrecht, 2002.
[5] Hurt, J.: Discriminant Analysis and Multidimensional Scaling. Wolfram Technology Conference. http//:library.wolfram.com/infocenter/conferences/5777. Champaing (IL), 2005.
[6] Statistical Methods in Finance. Handbook of Statistics. Volume 14. Maddala, G. S., Rao, C. R., eds. Elsevier. Amsterdam 1996.