﻿using System; namespace QuantRiskLib { ///Source: www.risk256.com /// ///References: ///Miller, Michael B. 2012. Mathematics and Statistics for Financial Risk Management. New York: John Wiley & Sons. ///Chapters 2, 3 and 10 public partial class Moments { #region Mean /// /// Returns the mean of an array. /// /// Array of data for which we are calculating the mean. public static double Mean(double[] array) { return Mean(array, 1.0); } /// /// Returns the mean of an array. /// /// Array of data for which we are calculating the mean. For time series, the last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in array is 1.0, the 2nd to last element d, 3rd to last d^2, ... public static double Mean(double[] array, double decayFactor) { int n = array.Length; double m = 0.0; double d = 1.0; for (int i = 0; i < n; i++) { m += d * array[n - 1 - i]; d *= decayFactor; } m *= GeometricSeries.InverseSumOfGeometricSeries(decayFactor, n); return m; } /// /// Returns the mean of an array. /// /// Array of data for which we are calculating the mean. For time series, the last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in array is 1.0, the 2nd to last element d, 3rd to last d^2, ... /// Window length. Method uses the most recent n points, n = length. public static double Mean(double[] array, double decayFactor, int length) { double[] subArray = Tools.MostRecentValues(array, length); return Mean(subArray, decayFactor); } /// /// Returns the weighted averages of the values in valueArray using the corresponding weights in weightArray. /// /// array of values for which we are computing the weighted average /// array of weights used in computing the weighted average public static double WeightedMean(double[] valueArray, double[] weightArray) { if (valueArray.Length != weightArray.Length) throw new ArgumentException("valueArray and weightArray must have the same number of elements"); double m = 0, s = 0; for (int i = 0; i < valueArray.Length; i++) { m += valueArray[i] * weightArray[i]; s += weightArray[i]; } m /= s; return m; } #endregion #region StandardDeviation /// /// Returns the sample standard deviation of an array. /// /// Array of data for which we are calculating the standard deviation. public static double StandardDeviation(double[] array) { return StandardDeviation(array, 1.0); } /// /// Returns the sample standard deviation of an array. /// /// Array of data for which we are calculating the standard deviation. For time series, the last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in array is 1.0, the 2nd to last element d, 3rd to last d^2, ... public static double StandardDeviation(double[] array, double decayFactor) { double v = Variance(array, decayFactor); return Math.Sqrt(v); } /// /// Returns the sample standard deviation of an array. /// /// Array of data for which we are calculating the standard deviation. For time series, the last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in array is 1.0, the 2nd to last element d, 3rd to last d^2, ... /// Window length. Method uses the most recent n points, n = length. public static double StandardDeviation(double[] array, double decayFactor, int length) { double[] subArray = Tools.MostRecentValues(array, length); return StandardDeviation(subArray, decayFactor); } /// /// Calculates the sample standard deviation of an array using a rolling window. /// The first element (index = 0) of the output array is the standard deviation of points  to [length - 1] in 'array', /// the second element of the output array is the standard deviation of points  to [length] in 'array', and so on. /// /// Array of data for which we are calculating the standard deviation. For time series, the last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in array is 1.0, the 2nd to last element d, 3rd to last d^2, ... /// >Length of rolling window length. Must be less than or equal to the lenght of the array of data. /// public static double[] RollingStandardDeviation(double[] array, double decayFactor, int length) { double[] rollingSdArray = RollingVariance(array, decayFactor, length); for (int i = 0; i < rollingSdArray.Length; i++) rollingSdArray[i] = Math.Sqrt(rollingSdArray[i]); return rollingSdArray; } #endregion #region Incremental Standard Deviation /// /// Returns the incremental standard deviation of a variable. Incremental standard deviation is to standard deviation what incremental VaR is to VaR. /// /// Return array of portfolio. The last element (index = n-1), is the most recent. /// Return array of the sub-portfolio for which incremental standard deviation is being measured. The last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in arrays is 1.0, the 2nd to last element d, 3rd to last d^2, ... public static double IncrementalStandardDeviation(double[] portfolioArray, double[] subPortfolioArray, double decayFactor) { if (portfolioArray.Length != subPortfolioArray.Length) throw new ArgumentException("Arrays must be the same length"); int len = portfolioArray.Length; double[] otherArray = new double[len]; for (int i = 0; i < len; i++) otherArray[i] = portfolioArray[i] - subPortfolioArray[i]; double varianceSubPort = Variance(subPortfolioArray, decayFactor); double covarSubOther = Covariance(subPortfolioArray, otherArray, decayFactor); //double sdPort = StandardDeviation(portfolioArray, decayFactor); double varianceOther = Variance(otherArray, decayFactor); double sdPort = Math.Sqrt(varianceSubPort + varianceOther + 2 * covarSubOther); double iStandardDeviation = (varianceSubPort + covarSubOther) / sdPort; return iStandardDeviation; } /// /// Returns the incremental standard deviation of a variable. Incremental standard deviation is to standard deviation what incremental VaR is to VaR. /// /// Return array of portfolio. The last element (index = n-1), is the most recent. /// Return array of the sub-portfolio for which incremental standard deviation is being measured. The last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in arrays is 1.0, the 2nd to last element d, 3rd to last d^2, ... /// Window length. Method uses the most recent n points, n = length. public static double IncrementalStandardDeviation(double[] portfolioArray, double[] subPortfolioArray, double decayFactor, int length) { double[] portfolioArrayLen = Tools.MostRecentValues(portfolioArray, length); double[] subPortfolioArrayLen = Tools.MostRecentValues(subPortfolioArray, length); return IncrementalStandardDeviation(portfolioArrayLen, subPortfolioArrayLen, decayFactor); } #endregion #region Variance /// /// Returns the sample variance of an array. /// /// Array of data for which we are calculating the variance. public static double Variance(double[] array) { return Variance(array, 1.0); } /// /// Returns the sample variance of an array. /// /// Array of data for which we are calculating the variance. For time series, the last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in array is 1.0, the 2nd to last element d, 3rd to last d^2, ... public static double Variance(double[] array, double decayFactor) { if(array.Length == 1) return 0.0; int n = array.Length; double v = 0.0; double d = 1.0; for (int i = 0; i < n; i++) { v += d * array[n - 1 - i] * array[n - 1 - i]; d *= decayFactor; } double m = Mean(array, decayFactor); double s1 = GeometricSeries.SumOfGeometricSeries(decayFactor, n); double s2 = GeometricSeries.SumOfGeometricSeries(decayFactor * decayFactor, n); double a = s1 / (s1 * s1 - s2); double b = s1 * a; v = a * v - b * m * m; //Occasionally floating-point error will cause v to be less than zero when it should be exactly zero or very close to zero. if(v < 0) { bool allEqual = Tools.ArrayAllEqual(array); return allEqual ? 0 : double.Epsilon; } return v; } /// /// Returns the sample variance of an array. /// /// Array of data for which we are calculating the variance. For time series, the last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in array is 1.0, the 2nd to last element d, 3rd to last d^2, ... /// Window length. Method uses the most recent n points, n = length. /// public static double Variance(double[] array, double decayFactor, int length) { double[] subArray = Tools.MostRecentValues(array, length); return Variance(subArray, decayFactor); } /// /// Calculates the sample variance of an array using a rolling window. /// The first element (index = 0) of the output array is the variance of points  to [length - 1] in 'array', /// the second element of the output array is the variance of points  to [length] in 'array', and so on. /// /// Array of data for which we are calculating the variance. For time series, the last element (index = n-1), is the most recent. /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in array is 1.0, the 2nd to last element d, 3rd to last d^2, ... /// >Length of rolling window length. Must be less than or equal to the lenght of the array of data. /// public static double[] RollingVariance(double[] array, double decayFactor, int length) { if (array.Length < length) throw new ArgumentException("Not enought points in array. Number of points in array must be greater than or equal to length."); double s1 = GeometricSeries.SumOfGeometricSeries(decayFactor, length); double s2 = GeometricSeries.SumOfGeometricSeries(decayFactor * decayFactor, length); double c = 1 / (s1 * s1 - s2); double a = s1 * c; double z1 = 0.0, z2 = 0.0; double d = 1.0; for (int i = 0; i < length; i++) { z1 += d * array[length - 1 - i] * array[length - 1 - i]; z2 += d * array[length - 1 - i]; d *= decayFactor; } double dToLength = Math.Pow(decayFactor, length); double[] vArray = new double[array.Length - length + 1]; vArray = a * z1 - c * z2 * z2; for(int i = 1; i < vArray.Length; i++) { z1 = decayFactor * z1 + array[length - 1 + i] * array[length - 1 + i] - dToLength * array[i - 1] * array[i - 1]; z2 = decayFactor * z2 + array[length - 1 + i] - dToLength * array[i - 1]; vArray[i] = a * z1 - c * z2 * z2; } return vArray; } #endregion #region Covariance and Correlation /// /// Returns the sample covariance between two arrays. /// Arrays should be of equal length, and contain more than one element. /// /// /// public static double Covariance(double[] array1, double[] array2) { if (array1.Length != array2.Length) throw new ArgumentException("Arrays must be the same length"); int n = array1.Length; if (n == 1) return double.NaN; double m1 = Mean(array1); double m2 = Mean(array2); double covar = 0.0; for (int i = 0; i < n; i++) covar += (array1[i] - m1) * (array2[i] - m2); covar /= (n - 1); return covar; } /// /// Returns the sample covariance between two arrays. /// Arrays should be of equal length, and contain more than one element. /// /// /// /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in arrays is 1.0, the 2nd to last element d, 3rd to last d^2, ... public static double Covariance(double[] array1, double[] array2, double decayFactor) { if (array1.Length != array2.Length) throw new ArgumentException("Arrays must be the same length"); if (Tools.ArrayAllEqual(array1) || Tools.ArrayAllEqual(array2)) return 0.0; int n = array1.Length; if (n == 1) return double.NaN; double covar = 0.0; double d = 1.0; for (int i = 0; i < n; i++) { covar += d * array1[n - 1 - i] * array2[n - 1 - i]; d *= decayFactor; } double m1 = Mean(array1, decayFactor); double m2 = Mean(array2, decayFactor); double S1 = GeometricSeries.SumOfGeometricSeries(decayFactor, n); double S2 = GeometricSeries.SumOfGeometricSeries(decayFactor * decayFactor, n); double A = S1 / (S1 * S1 - S2); double B = S1 * A; covar = A * covar - B * m1 * m2; return covar; } /// /// Returns the correlation between two arrays. /// Arrays should be of equal length, and contain more than one element. /// /// /// public static double Correlation(double[] array1, double[] array2) { double covar = Covariance(array1, array2); double sd1 = StandardDeviation(array1); double sd2 = StandardDeviation(array2); double corr = covar / (sd1 * sd2); //In theory we shouldn't get values less than -1 or greater than +1, but with precision errors we can get values slightly outside these bounds. if (corr < -1) return -1; if (corr > 1) return 1; return corr; } /// /// Returns the correlation between two arrays. /// Arrays should be of equal length, and contain more than one element. /// /// /// /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in arrays is 1.0, the 2nd to last element d, 3rd to last d^2, ... public static double Correlation(double[] array1, double[] array2, double decayFactor) { double covar = Covariance(array1, array2, decayFactor); double sd1 = StandardDeviation(array1, decayFactor); double sd2 = StandardDeviation(array2, decayFactor); double corr = covar / (sd1 * sd2); //In theory we shouldn't get values less than -1 or greater than +1, but with precision errors we can get values slightly outside these bounds. if (corr < -1) return -1; if (corr > 1) return 1; return corr; } /// /// Returns the serial covariance of an array of data. /// The serial covariance is the covariance of the data, with its own first lag. /// public static double SerialCovariance(double[] array) { int n = array.Length - 1; double[] sub1 = new double[n]; double[] sub2 = new double[n]; for (int i = 0; i < n; i++) { sub1[i] = array[i]; sub2[i] = array[i + 1]; } return Covariance(sub1, sub2); } /// /// Returns the serial covariance of an array of data. /// The serial covariance is the covariance of the data, with its own first lag. /// /// /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in arrays is 1.0, the 2nd to last element d, 3rd to last d^2, ... public static double SerialCovariance(double[] array, double decayFactor) { int n = array.Length - 1; double[] sub1 = new double[n]; double[] sub2 = new double[n]; for (int i = 0; i < n; i++) { sub1[i] = array[i]; sub2[i] = array[i + 1]; } return Covariance(sub1, sub2, decayFactor); } /// /// Returns the serial correlation of an array of data. /// The serial correlation is the correlation of the data, with its own first lag. /// public static double SerialCorrelation(double[] array) { int n = array.Length - 1; double[] sub1 = new double[n]; double[] sub2 = new double[n]; for (int i = 0; i < n; i++) { sub1[i] = array[i]; sub2[i] = array[i + 1]; } return Correlation(sub1, sub2); } /// /// Returns the serial correlation of an array of data. /// The serial correlation is the correlation of the data, with its own first lag. /// /// /// In most applications, the decay factor is between 0 and 1. Weigth on the last element in arrays is 1.0, the 2nd to last element d, 3rd to last d^2, ... public static double SerialCorrelation(double[] array, double decayFactor) { int n = array.Length - 1; double[] sub1 = new double[n]; double[] sub2 = new double[n]; for (int i = 0; i < n; i++) { sub1[i] = array[i]; sub2[i] = array[i + 1]; } return Correlation(sub1, sub2, decayFactor); } #endregion #region Skewness /// /// Returns the sample skewness of an array. /// /// Array of data for which we are calculating the skewness. public static double Skewness(double[] array) { int n = array.Length; if(n < 3) return double.NaN; double m = Mean(array); double sd = StandardDeviation(array); double skew = 0.0; for (int i = 0; i < n; i++) skew += Math.Pow(array[i] - m, 3.0); skew /= Math.Pow(sd, 3.0); skew = n * skew / ((n - 1.0) * (n - 2.0)); return skew; } #endregion #region Kurtosis /// /// Returns the sample kurtosis of an array. /// Note: This is not excess kurtosis. /// /// Array of data for which we are calculating the kurtosis. public static double Kurtosis(double[] array) { int n = array.Length; if (n < 4) return double.NaN; double m = Mean(array); double variance = Variance(array); double kurt = 0.0; for (int i = 0; i < n; i++) kurt += Math.Pow(array[i] - m, 4.0); kurt /= (variance * variance); kurt *= (n * (n + 1.0)) / ((n - 1.0) * (n - 2.0) * (n - 3.0)); return kurt; } /// /// Returns the sample excess kurtosis of an array. /// /// Array of data for which we are calculating the excess kurtosis. public static double ExcessKurtosis(double[] array) { double kurt = Kurtosis(array); int n = array.Length; kurt -= 3 * (n - 1.0) * (n - 1.0) / ((n - 2.0) * (n - 3.0)); return kurt; } #endregion #region Coskewness and Cokurtosis public enum CoskewnessType { AAB, ABB } public enum CokurtosisType { AAAB, AABB, ABBB } /// /// Returns the sample coskewness between two arrays. /// Arrays should be of equal length, and contain more than one element. /// /// /// /// public static double Coskewness(double[] arrayA, double[] arrayB, CoskewnessType coskewnessType) { double m3 = ThirdCentralCrossMoment(arrayA, arrayB, coskewnessType); double sd1 = StandardDeviation(arrayA); double sd2 = StandardDeviation(arrayB); if (coskewnessType == CoskewnessType.AAB) return m3/(sd1*sd1*sd2); return m3 / (sd1 * sd2 * sd2); } /// /// Returns the sample cokurtosis between two arrays. /// Arrays should be of equal length, and contain more than one element. /// /// /// /// public static double Coskurtosis(double[] arrayA, double[] arrayB, CokurtosisType cokurtosisType) { double m4 = FourthCentralCrossMoment(arrayA, arrayB, cokurtosisType); double sd1 = StandardDeviation(arrayA); double sd2 = StandardDeviation(arrayB); if (cokurtosisType == CokurtosisType.AAAB) return m4 / (sd1 * sd1 * sd1 * sd2); if(cokurtosisType == CokurtosisType.AABB) return m4 / (sd1 * sd1 * sd2 * sd2); return m4 / (sd1 * sd2 * sd2 * sd2); } /// /// Returns the sample third central cross moment between two arrays. /// Arrays should be of equal length, and contain more than one element. /// /// /// /// private static double ThirdCentralCrossMoment(double[] arrayA, double[] arrayB, CoskewnessType coskewnessType) { if(arrayA.Length != arrayB.Length) throw new ArgumentException("arrayA and arrayB should be the same length."); int n = arrayA.Length; double meanA = Mean(arrayA); double meanB = Mean(arrayB); double m3 = 0; if (coskewnessType == CoskewnessType.AAB) { for (int i = 0; i < n; i++) m3 += (arrayA[i] - meanA) * (arrayA[i] - meanA) * (arrayB[i] - meanB); } else { for (int i = 0; i < n; i++) m3 += (arrayA[i] - meanA) * (arrayB[i] - meanB) * (arrayB[i] - meanB); } m3 *= n / ((n - 1.0) * (n - 2.0)); //because this is the sample coskew multiply by n/(n-1)(n-2) rather than 1/n return m3; } /// /// Returns the sample fourth central cross moment between two arrays. /// Arrays should be of equal length, and contain more than one element. /// /// /// /// private static double FourthCentralCrossMoment(double[] arrayA, double[] arrayB, CokurtosisType cokurtosisType) { if (arrayA.Length != arrayB.Length) throw new ArgumentException("arrayA and arrayB should be the same length."); int n = arrayA.Length; double meanA = Mean(arrayA); double meanB = Mean(arrayB); double m4 = 0; if (cokurtosisType == CokurtosisType.AAAB) { for (int i = 0; i < n; i++) m4 += (arrayA[i] - meanA) * (arrayA[i] - meanA) * (arrayA[i] - meanA) * (arrayB[i] - meanB); } else if (cokurtosisType == CokurtosisType.AABB) { for (int i = 0; i < n; i++) m4 += (arrayA[i] - meanA)*(arrayA[i] - meanA)*(arrayB[i] - meanB)*(arrayB[i] - meanB); } else { for (int i = 0; i < n; i++) m4 += (arrayA[i] - meanA) * (arrayB[i] - meanB) * (arrayB[i] - meanB) * (arrayB[i] - meanB); } m4 *= n * (n - 1) / ((n - 1.0) * (n - 2.0) * (n - 3.0)); //because this is the sample coskew multiply by n(n+1)/(n-1)(n-2)(n-3) rather than 1/n return m4; } #endregion } } //Disclaimer //This code is freeware. The methods are not proprietary. Feel free to use, modify and redistribute. That said, if you plan //to use or redistribute give credit where credit is due and provide a link back to Risk256.com (or don't remove the link //and references already in the code). The code is intended primarily as an educational tool. No warranty is made as to the //code's accuracy. Use at your own risk.