Probability and Statistical Analysis

The mean absolute deviation (MAD) serves as a reliable measure of dispersion, reflecting the spread or variability of data points around the mean. However, when used in statistical inference, its statistical properties are inferior to the standard deviation. The mad function by default calculates the mean absolute deviation, and the formula is as follows:

The following results show that the calculations from DolphinDB and Python Numpy are consistent up to 8 decimal places.

import time
import numpy as np
import dolphindb as ddb

def calculate_mad(data):
    med = np.median(data)
    abs_diffs = np.abs(data - med)
    return np.median(abs_diffs)

# Set random seed to ensure reproducibility
np.random.seed(42)

data = np.random.randint(1, 101, size=10000)

# Calculate MAD using Python Numpy
result1 = calculate_mad(data)
print(f"MAD value using Python Numpy: {result1:.8f}")

# Calculate MAD using DolphinDB
s = ddb.session()
s.connect("localhost", 7271)
result2 = s.run("mad{,true}", data)
print(f"MAD value using DolphinDB: {result2:.8f}")

# Data validation
abs_diff = abs(result1 - result2)
if abs_diff < 1e-8:
    print("MAD calculation results are consistent up to 8 decimal places.")
else:
    print(f"MAD calculation results differ by: {abs_diff:.8f}")
    
s.close()

/* output
MAD value using Python Numpy: 25.00000000
MAD value using DolphinDB: 25.00000000
MAD calculation results are consistent up to 8 decimal places.
*/

The median absolute deviation (MAD) is a robust measure of statistical dispersion, providing an alternative to the mean absolute deviation that is less influenced by extreme outliers. It improves upon the traditional absolute deviation by using the median as a reference. To calculate the median absolute deviation, the mad function must specify useMedian=true. The formula for calculating MAD is as follows:

The following results show that the calculations from DolphinDB and Python Numpy are consistent up to 8 decimal places.

import time
import numpy as np
import dolphindb as ddb

def calculate_mad(data):
    med = np.median(data)
    abs_diffs = np.abs(data - med)
    return np.median(abs_diffs)

# Set random seed to ensure reproducibility
np.random.seed(42)

data = np.random.randint(1, 101, size=10000)

# Calculate MAD using Python Numpy
result1 = calculate_mad(data)
print(f"MAD value using Python Numpy: {result1:.8f}")

# Calculate MAD using DolphinDB
s = ddb.session()
s.connect("localhost", 7271)
result2 = s.run("mad{,true}", data)
print(f"MAD value using DolphinDB: {result2:.8f}")

# Data validation
abs_diff = abs(result1 - result2)
if abs_diff < 1e-8:
    print("MAD calculation results are consistent up to 8 decimal places.")
else:
    print(f"MAD calculation results differ by: {abs_diff:.8f}")
    
s.close()

/* output
MAD value using Python Numpy: 25.00000000
MAD value using DolphinDB: 25.00000000
MAD calculation results are consistent up to 8 decimal places.
*/

MAD can be effectively used to detect outliers. The process typically involves the following steps:

1. Use the med function to obtain the median of the dataset to serve as the central reference.
2. For each data point, use the mad function to calculate the median absolute deviation.
3. Set an outlier threshold, usually as a multiple of MAD.

   

4. Detect and mark outliers for further analysis.

   

The following example demonstrates how to use median absolute deviation for identifying and handling outliers via winsorization:

/* Identifying outliers based on MAD
   Input : original data
   Output : indexes of outliers
 */
def winsorized(X){
    // calculate the median and mean absolute deviation
    medianX = median(X)
    madX = mad(X, true)
    sigma = 1.4826*madX
    // calculate the lower and upper limits
    lowerLimit = medianX - 3*sigma
    upperLimit = medianX + 3*sigma
    // winsorize the data
    indexes = 0..(X.size()-1)
    return indexes.at(X < lowerLimit || X > upperLimit)
}
X = [-5.8, 2, 3, 4, 5, 6, 7, 8, 9, 18]
winsorized(X)
// output: [0, 9]

Skewness measures the asymmetry of the probability distribution of a real-valued random variable about its mean. The skewness value can be:

• Positive (right-skewed): The tail is more pronounced on the right side than on the left.
• Negative (left-skewed): The tail is more pronounced on the left rather than the right.
• Zero: Values are relatively evenly distributed around the mean.

Skewness can be computed using the skew function when the sample size n≥3. Let x̄ be the sample mean, s be the sample standard deviation, and m₃ be the third central moment of the sample. The formula for sample skewness is:

Note that the skewness computed in Python and Excel differs from the above. These tools calculate the Fisher-Pearson standardized moment G₁, which adjusts the sample skewness for bias. The formula is:

To compute the Fisher-Pearson standardized moment G₁ in DolphinDB, you can use the following script. The results show that the calculations from DolphinDB and Python pandas are consistent up to 8 decimal places.

import pandas as pd
import numpy as np
import dolphindb as ddb

# Set random seed to ensure reproducibility
np.random.seed(42)
data = np.random.normal(0, 1, 1000)

# Calculate Fisher-Pearson skewness using Python pandas
result1 = pd.Series(data).skew()
print(f"skewness using Python pandas: {result1:.8f}")

# Calculate Fisher-Pearson skewness using DolphinDB
session = ddb.session()
session.connect("localhost", 7271)
session.upload({"n":len(data),"x":data})
result2 = session.run("sqrt(n*(n-1))\(n-2)*skew(x)")
print(f"skewness using DolphinDB: {result2:.8f}")

# Data validation
abs_diff = abs(result1 - result2)
if abs_diff < 1e-8:
    print("Skewness calculation results are consistent up to 8 decimal places.")
else:
    print(f"Skewness calculation results differ by: {abs_diff:.8f}")

session.close()

/* output
Skewness value using Python pandas: 0.11697637
Skewness value using DolphinDB: 0.11697637
Skewness calculation results are consistent up to 8 decimal places.
*/

Kurtosis refers to a descriptive statistic used to measure how data disperse between the center and tails of a distribution. A higher kurtosis indicates that the variance is primarily driven by extreme values (either larger or smaller than the mean) that occur less frequently.

In DolphinDB, the kurtosis is calculated as an estimate of the fourth central moment. Subtracting 3 from the result yields a value that aligns with Python’s kurtosis calculation.

• By default, biased=true is used to calculate the sample kurtosis, which results in a biased estimate. Let the sample mean be denoted by x̄, the sample standard deviation by s, and the fourth central moment by m₄. The formula for sample kurtosis is as follows:

  

• If biased=false is specified, the calculation is for the population kurtosis, which provides an unbiased estimate for a normal distribution. The formula is as follows:

  

Since the kurt function in Python pandas calculates the unique symmetric unbiased estimator of the fourth cumulant (i.e., excess kurtosis), the results from DolphinDB's kurtosis function match those of pandas after subtracting 3, as shown in the following script.

import pandas as pd
import Numpy as np
import dolphindb as ddb

# Set random seed to ensure reproducibility
np.random.seed(42)
data = np.random.normal(0, 1, 1000)

# Calculate kurtosis using Python pandas
result1 = pd.Series(data).kurt()
print(f"Kurtosis (Pandas): {kurt_pandas:.8f}")

# Calculate kurtosis using DolphinDB
session = ddb.session()
session.connect("localhost", 7271)
result2 = session.run("kurtosis{,false}",data)-3
print(f"Kurtosis (DolphinDB): {kurt_dolphindb:.8f}")

# Data validation
abs_diff = abs(result1 - result2)
if abs_diff < 1e-8:
    print("Kurtosis calculation results are consistent up to 8 decimal places.")
else:
    print(f"Kurtosis calculation results differ by: {abs_diff:.8f}")

session.close()

/* output
Kurtosis (Pandas): 0.07256220
Adjusted Kurtosis (DolphinDB): 0.07256220
Kurtosis calculation results are consistent up to 8 decimal places.
*/

Covariance measures the joint variability of two random variables. A covariance matrix is a square matrix giving the covariance between each pair of elements across all dimensions. In DolphinDB, the covar function calculates the sample covariance between two vectors using the following formula:

For multi-dimensional data, DolphinDB provides the covarMatrix function to compute the covariance matrix, using the formula:

The following example demonstrates how to construct a risk matrix in DolphinDB:

returns_a = [0.01, 0.02, -0.03, 0.01, 0.02]
returns_b = [-0.02, 0.01, 0.03, -0.02, -0.01]
returns_c = [0.03, -0.01, -0.02, 0.02, 0.01]
m = matrix(returns_a, returns_b, returns_c).rename!(`a`b`c)
covariance_matrix = covarMatrix(m)
/* output:
	a	b	c
a	0.0004300000	-0.0003100000	0.0002300000
b	-0.0003100000	0.0004700000	-0.0004350000
c	0.0002300000	-0.0004350000	0.0004300000
*/

Pearson correlation coefficient measures linear correlation between two variables. In DolphinDB, the correlation coefficient can be calculated using the corr function, using the formula:

For multi-dimensional data, DolphinDB provides the corrMatrix function to compute the correlation matrix, using the formula:

The Pearson correlation coefficient between a market index and individual stock returns shows the correlation between the stock and the overall market. For example:

// Market index data
market_index = [100, 110, 120, 130, 140]
// Individual stock return data
stock_returns = [0.01, 0.02, 0.03, -0.01, 0.02]
corr(market_index, stock_returns)
// Output: -0.1042572070

In portfolio optimization, correlation matrix can help reveal the correlation between different assets. This allows investors to select assets with low correlations in the portfolio for better diversification. For example, select assets with low correlations to others:

setRandomSeed(42)
m = rand(10.0, 100)$10:10
c= corrMatrix(m)
result = table(1..10 as stock_code, 
               c[c<0.5].abs().sum() as corr_sum).sortBy!(`corr_sum)
result
/* Output
	stock_code	corr_sum
0	8	1.6943835693
1	2	1.9124975741
2	6	1.9699810000
3	10	1.9706433980
4	4	1.9955912515
5	1	2.0791768735
6	9	2.2377748096
7	5	2.2757034717
8	7	2.8368851661
9	3	3.0716024528
*/

The Spearman’s rank correlation coefficient is a non-parametric measure of rank correlation, essentially the Pearson correlation coefficient applied to ranked variables. In DolphinDB, it can be computed using the spearmanr function, using the formula:

In quantitative trading, the Information Coefficient (IC) measures the correlation between predicted and actual future returns, often used to evaluate the forecasting and stock-picking abilities. The spearmanr function is commonly used to calculate the Rank IC values. The result of spearmanr is consistent with the rank-based correlation computed using corr.

setRandomSeed(123)
predicted_returns = rand(1.0, 100)
actual_returns = rand(1.0, 100)
r1 = spearmanr(predicted_returns, actual_returns)
r2 = corr(rank(predicted_returns) + 1, rank(actual_returns) + 1)
eqObj(r1, r2)
// output: true

Bernoulli distribution (0-1 distribution) is a basic univariate distribution and serves as the basis for the Binomial distribution. Suppose an investment strategy generates daily returns as a binary random variable, which could either be a positive return (1) or a negative return (0), with no cases of no gain or loss. This aligns with the following description:

The probability of positive returns p_win determines all the characteristics of the Bernoulli distribution. For example, if p_win=0.5, the probability mass function (PMF) of the 0-1 distribution can be calculated using the pmfBernoulli function in DolphinDB.

def pmfBernoulli(p){
   cdf = cdfBinomial(1, p, [0,1])
   pmf =table([0,1] as X,cdf -prev(cdf).nullFill(0) as P).sortBy!(`X)
   return pmf
}
p = 0.5
pmfBernoulli(p)
/* Output 
	X	P
0	0	0.50
1	1	0.50
*/

In DolphinDB, once the probability p_win of a 0-1 event and the number of trials N are defined, the randBernoulli function can be used to simulate Bernoulli trials. For instance, assume each positive return corresponds to a logarithmic return of 10%, and each negative return corresponds to -10%. After running N trials according to the strategy, the total return can be calculated by summing the logarithmic returns. The following example shows that the return after 10 investment periods is 22%.

setRandomSeed(40393)
def randBernoulli(p, n){
   return randBinomial(1, p, n)
}
x = randBernoulli(0.5, 10)
// Suppose each period's logarithmic return is fixed at 0.1
log_return = 0.1
exp(sum(log_return * iif(x == 0, -1, x))) - 1
// Output: 0.2214

Evaluating the probability and frequency of positive returns can be modeled using the Binomial distribution. A random variable X that follows a Binomial distribution with parameters p (probability of a 0-1 event) and n (number of trials) is denoted as X∈B (n, p) . The probability mass function (PMF) for X=k is expressed as:

In DolphinDB, the pmfBinomial function can be used to compute the PMF of a Binomial distribution. For example, the following result shows that after 10 trials, the probability of achieving exactly 5 positive returns is approximately 24.61%.

def pmfBinomial(trials, p){
    X = 0..trials
    cdf = cdfBinomial(trials, p, X)
    return table(X, cdf - prev(cdf).nullFill(0) as P).sortBy!(`X)
}
pmfBinomial(10, 0.5)

/* output:
X  P                
-- -----------------
0  0.0009765625     
1  0.009765625      
2  0.0439453125     
3  0.1171875        
4  0.205078125      
5  0.246093749999999
6  0.205078125      
7  0.1171875        
8  0.0439453125     
9  0.009765625      
10 0.0009765625   
*/

Using the pmfBinomial function, we can visualize the Binomial distribution under varying probabilities p and trial numbers n. Note that only integer values of X are valid.

plot([pmfBinomial(20,0.5)["P"] <- take(double(),20) as "p=0.5,n=20",
      pmfBinomial(20,0.7)["P"] <- take(double(),20) as "p=0.7,n=20",
      pmfBinomial(40,0.5)["P"] as "p=0.5,n=40"],0..40,"Binomial Distributions")

The Poisson distribution is a statistical model used to describe the occurrence of events over a continuous space or time interval. It is closely related to the Binomial distribution but differs in focus. While the Binomial distribution models the number of successes in a fixed number of trials, the Poisson distribution focuses on the number of discrete events occurring within a continuous interval. The probability mass function (PMF) of the Poisson distribution is defined as:

In DolphinDB, the pmfPoisson function can be used to compute the PMF values of a Poisson distribution. Note that only integer values of X are valid, and visual grouping is often done with lines connecting discrete points.

def pmfPoisson(mean, upper=20){
   x = 0..upper
   cdf =cdfPoisson(mean, x)
   pmf =table(x as X, cdf - prev(cdf).nullFill(0) as P).sortBy!(`X)
   return pmf
}
plot([pmfPoisson(1)["P"] as "λ=1",
      pmfPoisson(4)["P"] as "λ=4", 
      pmfPoisson(10)["P"] as "λ=10"],
      0..20,"Poisson Distributions")

The t-distribution is widely used in statistics, particularly when making inferences when the population mean and variance is unknown. If x̄ is the sample mean and s is the sample standard deviation, the resulting t statistic is expressed as:

A common application of the t-distribution is in calculating the confidence interval for the mean. In DolphinDB, the critical value corresponding to the t statistic can be obtained using the invStudent function. For example, assuming a sample size of n=20, the degrees of freedom df=19, and a significance level of 95%, the critical value corresponding to t_α/2 can be calculated using invStudent(19, 0.975).

x = randNormal(-100, 1000, 10000)
n = x.size()
df = n-1
rt = invStudent(df, 0.975)
// output: 1.9602012636

If a random variable X∈N(μ , σ), then . The sum of squares of n independent standard normal random variables follows a Chi-Square distribution with n degrees of freedom:

When the sample size is small, the ratio of the sample standard deviation S² to the population standard deviation σ² also follows a Chi-Square distribution, expressed as:

For example, consider stock returns with a restricted volatility σ²=0.05. To analyze whether the standard deviation exceeds the permissible value, the critical value of the Chi-Square statistic can be calculated using DolphinDB's invChiSquare function. The analysis shows that the Chi-Square statistic does not exceed the 95% significance level critical value, so the null hypothesis S² <=σ² cannot be rejected. This implies that the sample return volatility is preliminarily inferred to be within the specified range.

setRandomSeed(123)
x = randUniform(-0.1, 0.1, 100)
test = (x.size() - 1) * pow(x.std(), 2) \ 0.05
test > invChiSquare(x.size() - 1, 0.95)
// output: false

F-Distribution

The F-distribution is named after Ronald Fisher, who introduced it to determine critical values in ANOVA (Analysis of Variance). By calculating the ratio of variances from two investment portfolios, the F-distribution helps assess whether the two portfolios exhibit the same risk level:

In this context, S_x² and S_y² represent the sample standard deviations of the stock returns from the two portfolios, and the distribution of their ratio is the F-distribution. In ANOVA applications, the critical value of the F-distribution is determined based on the degrees of freedom of the numerator and denominator and the significance level, as shown below:

In DolphinDB, the critical value of the F-distribution can be calculated using the invF function. Assuming that the two investment portfolios come from different Gaussian distributions, with the latter having a higher standard deviation, the compareTwoSamples function is used to perform an F-test. The result indicates that we can reject the null hypothesis σ₁ = σ₂, which matches the expected outcome.

setRandomSeed(1)
x1 = randNormal(0, 1, 20)
x2 = randNormal(0, 5, 30)

def compareTwoSamples(x1, x2, alpha){
   test = (pow(x1.std(), 2)\(x1.size()-1))\(pow(x2.std(), 2)\(x2.size()-1))
   return test > invF(x1.size()-1, x2.size()-1, 1-alpha\2) || 
          test < invF(x1.size()-1, x2.size()-1, alpha\2)
}

compareTwoSamples(x1, x2, 0.05)
// output: true

Using the data from Chapter 4 as an example, the time series graph of stock returns is plotted below.

setRandomSeed(123)
facTable = getAllFactorTable(2, 240)
tmp = select return_day from facTable pivot by record_date, stock_code
plot(tmp[,1:], tmp["record_date"])

We observe that the returns for both stocks are nearly identical, suggesting that their mean returns may be the same. Next, we perform a hypothesis test to verify this assumption:

Hypothesis: The return samples for both stocks follow a normal distribution, with known variances, and both stocks have a standard deviation of 6.

We use the zTest function to perform the z-test, and the output is as follows. The calculated z-statistic's p-value is 0.68 (> 0.05), meaning we fail to reject the null hypothesis, and conclude that the mean returns of the two stocks are indeed the same.

zTest(tmp["600001"], tmp["600002"], 0, 6, 6)
/* output:
method: 'Two sample z-test'
zValue: 0.4050744166
confLevel: 0.9500000000
stat:
alternativeHypothesis	pValue	lowerBound	upperBound
0	difference of mean is not equal to 0	0.6854233525	-0.8516480908	1.2953848817
1	difference of mean is less than 0	0.6572883237	-∞	1.1227918307
2	difference of mean is greater than 0	0.3427116763	-0.6790550398	∞
*/

In practice, since the volatility of stock returns is unknown and the sample size is typically small, the t-test should be used instead of z-test. Using DolphinDB's tTest function, the output is as follows. The calculated t-statistic's p-value is 0.13 (> 0.05), so we fail to reject the null hypothesis, concluding that the mean returns of the two stocks are indeed the same.

tTest(tmp["600001"], tmp["600002"], 0)
/* output:
method: 'Welch two sample t-test'
tValue: -0.6868291027
df: 473.7980010389
confLevel: 0.9500000000
stat:
    alternativeHypothesis	pValue	lowerBound	upperBound
0	difference of mean is not equal to 0	0.1314	-2.0292	0.2651
1	difference of mean is less than 0	0.0657	-∞	0.0801
2	difference of mean is greater than 0	0.9342	-1.8441	∞
*/

Using the dataset mentioned above, test if the volatility (variance) of stock 600001's returns is greater than a specific threshold (6). The null and alternative hypotheses are:

The chi-square statistic is calculated using the following formula, and the critical value at the 95% significance level is obtained using DolphinDB's invChiSquare function. The results show that the chi-square statistic is less than the critical value, so we fail to reject the null hypothesis, concluding that the volatility of stock 600001 is not greater than 6.

test_statistic =(size(tmp)-1)*tmp["600001"].var()/36
print(stringFormat("chi-squre test statistic:%.2F",test_statistic))
chi_statistic=invChiSquare(size(tmp)-1,0.95)
print(stringFormat("chi-squre test statistic:%.2F",chi_statistic))
/* output:
chi-squre test statistic:233.59
chi-squre test statistic:276.06
*/

To compare whether the volatilities (variances) of returns for two stocks are consistent, an F-test can be used. The null and alternative hypotheses are as follows:

Using DolphinDB's fTest function, the F-test output is as follows. Based on the analysis, the calculated F-statistic is 0.79, which lies between the left and right critical values of the F-distribution. Therefore, we fail to reject the null hypothesis, concluding that the variances of the returns for the two stocks are equal.

fTest(tmp["600001"], tmp["600002"])
/*Output 
method: 'F test to compare two variances'
fValue: 0.7907099827
numeratorDf: 239
denominatorDf: 239
confLevel: 0.9500000000
stat:
	alternativeHypothesis	pValue	lowerBound	upperBound
0	ratio of variances is not equal to 1	0.0701139943	0.6132461206	1.0195291184
1	ratio of variances is less than 1	0.0350569972	0.0000000000	0.9786251182
2	ratio of variances is greater than 1	0.9649430028	0.6388782232	∞
*/