/*****************************************************************************
*                                                                            *
*   KURTOSIS.SPL Copyright (C) 2008 DSP Development Corporation              *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Calculates the kurtosis of a distribution                   *
*                                                                            *
*   Revisions:    6 Feb 2008  RRR  Creation                                  *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_KURTOSIS

    KURTOSIS

    Purpose: Calculates the kurtosis of a series.

    Syntax:  KURTOSIS(s, biased, excess, vardef)

                     s - Optional series. Defaults to the current
                         Window.

                biased - Optional, an integer.  Bias calculation,
                           0: compute unbiased population estimate (default)
                           1: biased estimate

                excess - Optional, a real. Excess normalization subtracted
                         from the base calculation. Defaults to 3.0.


                vardef - Optional, an integer. Variance normalization method.
                           0: 1/(N-1)
                           1: 1/N

                           where N is the length of the input series.
                           Defaults to 1/N for the biased computation
                           and 1/(N-1) for the unbiased computation.

    Returns: A scalar, the kurtosis of the series

    Example:
             W1: {3, 4, 5, 2, 3, 4, 5, 6, 4, 7}
             W2: {kurtosis(w1)}
             W3: {kurtosis(w1, 1)}

             W2 == {-0.151800}, the unbiased or population kurtosis.
             W3 == {-0.631321}, the biased kurtosis.

    Remarks:
             Kurtosis characterizes the relative peakedness or flatness
             of a distribution compared with the normal distribution.
             Positive kurtosis indicates a relatively peaked
             distribution. Negative kurtosis indicates a relatively
             flat distribution.

             A distribution with a high peak is called leptokurtic, a
             flat-topped curve is called platykurtic, and the normal
             distribution is called mesokurtic.

             By default, KURTOSIS computes the unbiased or population
             kurtosis as used in SAS, SPSS, and Excel. To calculate the
             biased kurtosis, set the BIASED flag to 1.

             By default, KURTOSIS subtracts 3.0 from the raw kurtosis
             computation since the kurtosis of a normal distribution is 3.
             Set EXCESS to 0.0 to compute the raw kurtosis.

             Matlab computes the raw biased kurtosis by default. To match
             the results of Matlab for either biased or unbiased kurtosis,
             add 3.0. For example:

                ml.kurtosis(s)    == kurtosis(s, 1) + 3;
                ml.kurtosis(s, 1) == kurtosis(s, 1) + 3;
                ml.kurtosos(s, 0) == kurtosis(s, 0) + 3;

             For multi-column data, KURTOSIS returns the kurtosis of each
             column.

    See Also:
             Skew
             Stdev
#endif


/* calculate biased or unbiased kurtosis of a series */
ITERATE kurtosis(s, biased, excess, vardef)
{
        local k, n;

        (s, biased, excess, vardef) = kurtosis_parse_args(s, biased, excess, vardef);

        n = length(s);

        if (vardef == -1)
        {
                vardef = (biased) ? 1 : 0;
        }

        /* base calculation */
        k = sum(abs(s - mean(s)) ^ 4) / std(s, vardef) ^ 4;

        if (biased)
        {
                k = k / n - excess;
        }
        else
        {
                /* unbiased kurtosis population */
                if (n > 3)
                {
                        k *= (n * (n + 1)) / ((n - 1) * (n - 2) * (n - 3));
                        k -= (excess * (n - 1) ^ 2) / ((n - 2) * (n - 3));
                }
                else
                {
                        k = nan;
                }
        }
        
        return(k);
}


kurtosis_parse_args(s, biased, excess, vardef)
{
        local k, n;

        if (argc < 4)
        {
                if (argc < 3)
                {
                        if (argc < 2)
                        {
                                if (argc < 1)
                                {
                                        s      = refseries(w0);
                                        biased = 0;
                                }
                                else if (not(isarray(s)))
                                {
                                        biased = s;
                                        s      = refseries(w0);
                                }
                                else
                                {
                                        biased = 0;
                                }
                        }
                        if (not(isarray(s)))
                        {
                                excess = biased;
                                biased = s;
                                s      = refseries(w0);
                        }
                        else
                        {
                                excess = 3.0;
                        }
                }
                
                vardef = -1;
        }
        
        return(s, biased, excess, vardef);
}