For financial data I have successfully used heavy-tail Lambert W x Gaussian transformations.
Pyhon: gaussianize
is an sklearn-type implementation of the IGMM
algorithm in Python.
C++: the lamW R package has an elegant (and fast) C++ implementation of Lambert's W function. This can be a starting point for a full C++ implementation of IGMM or MLE for Lambert W x Gaussian transformations.
R: the LambertW R package is a full implementation of the Lambert W x F framework (simulation, estimation, plotting, transformation, testing).
As an illustration consider the SP500
return series in R.
library(MASS)
data(SP500)
yy <- ts(SP500)
library(LambertW)
test_norm(yy)
## $seed
## [1] 516797
##
## $shapiro.wilk
##
## Shapiro-Wilk normality test
##
## data: data.test
## W = 1, p-value <2e-16
##
##
## $shapiro.francia
##
## Shapiro-Francia normality test
##
## data: data.test
## W = 1, p-value <2e-16
##
##
## $anderson.darling
##
## Anderson-Darling normality test
##
## data: data
## A = 20, p-value <2e-16
As is well-known, financial data typically have fat tails and are sometimes negatively skewed. For the SP500
case, skewness is not too large, but it exhibits high kurtosis (7.7). Also several normality tests clearly reject the null hypothesis of a marginal Gaussian distribution.
Since we only have to deal with heavy tails, but not skewness, let's fit a heavy-tailed Lambert W x Gaussian distribution using a method of moments estimator (one could also use the maximum likelihood estimator (MLE) with MLE_LambertW()
).
# fit a heavy tailed Lambert W x Gaussian
mod <- IGMM(yy, type = "h")
mod
## Call: IGMM(y = yy, type = "h")
## Estimation method: IGMM
## Input distribution: Any distribution with finite mean & variance and kurtosis = 3.
## mean-variance Lambert W x F type ('h' same tails; 'hh' different tails; 's' skewed): h
##
## Parameter estimates:
## mu_x sigma_x delta
## 0.05 0.72 0.16
##
## Obtained after 4 iterations.
The heavy tail parameter $\widehat{\delta} = 0.16$ is significantly different from zero and implies heavy tails. For such a $\delta$ moments up to order $1 / \widehat{\delta} = 6.29$ exist.
The model check question is of course if the back-transformed data does indeed have a Gaussian distribution. Let's check again using test_norm()
:
# transform data to input data (which presumably should have Normal distribution); use return.u = TRUE to get zero-mean, unit variance data
xx <- get_input(mod, return.u = FALSE)
test_norm(xx)
## $seed
## [1] 268951
##
## $shapiro.wilk
##
## Shapiro-Wilk normality test
##
## data: data.test
## W = 1, p-value = 0.2
##
##
## $shapiro.francia
##
## Shapiro-Francia normality test
##
## data: data.test
## W = 1, p-value = 0.2
##
##
## $anderson.darling
##
## Anderson-Darling normality test
##
## data: data
## A = 0.7, p-value = 0.07
I think the plot and normality test results speak for themselves.
The package also provides a single function that does all these steps at once: Gaussianize()
(this is also what the Python package implements).