/*****************************************************************************
*                                                                            *
*   COVM.SPL     Copyright (C) 1998-2017 DSP Development Corporation         *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:      Randy Race                                                  *
*                                                                            *
*   Synopsis:    Calculates the covariance matrix                            *
*                                                                            *
*   Revisions:   19 Feb 1998  RRR  Creation - from COVM.MAC                  *
*                 3 May 2000  RRR  optional second array                     *
*                 2 Mar 2017  RRR  array shaping for 2 inputs                *
*                                                                            *
*****************************************************************************/

#if @HELP_COVM

    COVM

    Purpose: Calculates the covariance matrix of an array

    Syntax:  COVM(m1, m2)

              m1 - an array or series

              m2 - an optional second array or series

    Returns: An array

    Example:


            a = {{1.00, 3.00, 2.20},
                 {1.10, 4.00, 2.40},
                 {1.20, 5.00, 2.60}}

            b = covm(a)

            b == {{0.01, 0.10, 0.02},
                  {0.10, 1.00, 0.20},
                  {0.02, 0.20, 0.04}}

            c = diag(sqrt(b))
            d = colstdev(a)'

            c == d


    Remarks:

            The mean is removed from each column before the
            covariance is computed.  The standard deviations of each
            column can be calculated by:

                      diag(sqrt(covm(m)))

    See Also
             *^
             Colstdev
#endif


/* calculate covariance matrix */
covm(m1, m2)
{
        local cov;

        if (argc < 2)
        {
                if (argc < 1) error("covm - input series or array required");
        }
        else
        {
                if (not(all(size(m1) == size(m2))))
                {
                        error("covm - input sizes do not match");
                }

                m1 = ravel(unravel(m1), unravel(m2));
        }

        try
        {
                /* fast real VSL covariance */
                cov = covariance(m1);
        }
        catch
        {
                /* remove mean from each column in a matrix */
                cov = (crossprod(m1 - colmean(m1)) / (length(m1) - 1));
        }

        return(cov);
}