Skip to content

mfe.distributions

mfe.distributions

mfe.distributions — Fat-tailed and multivariate distributions.

skewt_logpdf Hansen (1994) Skew-t log-PDF with analytic score skewt_ppf Skew-t quantile function (VaR/ES) ged_logpdf Generalized Error Distribution log-PDF ged_ppf GED quantile function mvnorm_loglik Multivariate normal log-likelihood (time-varying Sigma_t) mahalanobis Mahalanobis distances under a covariance sequence standardize_mvn Extract standardized multivariate residuals

skewt_logpdf 

skewt_logpdf(x: FloatArray, nu: float, lam: float) -> FloatArray

Log-PDF of Hansen's Skew-t distribution.

Parameters:

Name Type Description Default
x FloatArray
required
nu float
required
lam skewness(-1, 1)
required
Source code in src/mfe/distributions/skewt.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def skewt_logpdf(
    x: FloatArray,
    nu: float,
    lam: float,
) -> FloatArray:
    """
    Log-PDF of Hansen's Skew-t distribution.

    Parameters
    ----------
    x   : (T,) standardized residuals
    nu  : degrees of freedom (>2)
    lam : skewness (-1, 1)
    """
    x = np.asarray(x, dtype=np.float64)
    a, b, c = _hansen_constants(nu, lam)

    mask = x < -a / b
    z = np.where(mask,
                 (b * x + a) / (1 - lam),
                 (b * x + a) / (1 + lam))

    log_kernel = -(nu + 1) / 2 * np.log(1 + z ** 2 / (nu - 2))
    log_pdf = np.log(b) + np.log(c) + log_kernel

    return log_pdf

skewt_score 

skewt_score(x: FloatArray, nu: float, lam: float) -> tuple[FloatArray, FloatArray]

Analytic score of log-skewt PDF with respect to (nu, lam).

Returns (d_log_f / d_nu, d_log_f / d_lam), each (T,).

These are used in the outer MLE loop over distribution parameters, avoiding the finite-difference approach from the MATLAB source.

Source code in src/mfe/distributions/skewt.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def skewt_score(
    x: FloatArray,
    nu: float,
    lam: float,
) -> tuple[FloatArray, FloatArray]:
    """
    Analytic score of log-skewt PDF with respect to (nu, lam).

    Returns (d_log_f / d_nu, d_log_f / d_lam), each (T,).

    These are used in the outer MLE loop over distribution parameters,
    avoiding the finite-difference approach from the MATLAB source.
    """
    from scipy.special import digamma

    x = np.asarray(x, dtype=np.float64)
    a, b, c = _hansen_constants(nu, lam)

    mask = x < -a / b
    sign = np.where(mask, 1 - lam, 1 + lam)
    z = (b * x + a) / sign

    u = z ** 2 / (nu - 2)
    denom = 1 + u

    # d/d_nu: from log kernel
    d_kernel_dnu = (
        -0.5 * np.log(denom)
        + (nu + 1) / 2 * z ** 2 / ((nu - 2) ** 2 * denom)
    )
    d_c_dnu = 0.5 * (digamma((nu + 1) / 2) - digamma(nu / 2) - 1 / (nu - 2))
    score_nu = d_c_dnu + d_kernel_dnu

    # d/d_lam: via chain rule through z(lam) and a(lam), b(lam)
    # This is the expensive part — see Hansen (1994) Appendix
    da_dlam = 4 * c * (nu - 2) / (nu - 1)  # simplified: ignoring dc/dlam for now
    db_dlam = (6 * lam - 2 * a * da_dlam) / (2 * b)

    dz_dlam_pos = (db_dlam * x + da_dlam) * (1 + lam) - (b * x + a)
    dz_dlam_pos /= (1 + lam) ** 2
    dz_dlam_neg = (db_dlam * x + da_dlam) * (1 - lam) + (b * x + a)
    dz_dlam_neg /= (1 - lam) ** 2

    dz_dlam = np.where(mask, dz_dlam_neg, dz_dlam_pos)
    d_kernel_dlam = -(nu + 1) * z * dz_dlam / ((nu - 2) * denom)
    score_lam = db_dlam / b + d_kernel_dlam

    return score_nu, score_lam

skewt_ppf 

skewt_ppf(p: FloatArray, nu: float, lam: float) -> FloatArray

Quantile function (inverse CDF) of Hansen's Skew-t. Used for VaR/ES computation.

Source code in src/mfe/distributions/skewt.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def skewt_ppf(
    p: FloatArray,
    nu: float,
    lam: float,
) -> FloatArray:
    """
    Quantile function (inverse CDF) of Hansen's Skew-t.
    Used for VaR/ES computation.
    """
    p = np.asarray(p, dtype=np.float64)
    a, b, c = _hansen_constants(nu, lam)

    p1 = (1 - lam) / 2  # mass in the left tail

    result = np.where(
        p < p1,
        (student_t.ppf(p / (1 - lam), df=nu) * (1 - lam) - a) / b,
        (student_t.ppf((p - (1 - lam) / 2) / (1 + lam) + 0.5, df=nu) * (1 + lam) - a) / b,
    )
    return result

ged_logpdf 

ged_logpdf(x: FloatArray, nu: float) -> FloatArray

Log-PDF of the GED with zero mean and unit variance.

Parameters:

Name Type Description Default
x FloatArray
required
nu shape parameter (nu > 0); nu=2 is Normal
required
Source code in src/mfe/distributions/ged.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
def ged_logpdf(x: FloatArray, nu: float) -> FloatArray:
    """
    Log-PDF of the GED with zero mean and unit variance.

    Parameters
    ----------
    x  : (T,) standardized residuals
    nu : shape parameter (nu > 0); nu=2 is Normal
    """
    x = np.asarray(x, dtype=np.float64)
    lam = _ged_lambda(nu)

    log_c = np.log(nu) - np.log(2 * lam) - np.log(gamma(1 / nu))
    log_kernel = -0.5 * np.abs(x / lam) ** nu

    return log_c + log_kernel

ged_ppf 

ged_ppf(p: FloatArray, nu: float) -> FloatArray

Quantile function of the GED(nu). Uses the relationship to the gamma distribution.

Source code in src/mfe/distributions/ged.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def ged_ppf(p: FloatArray, nu: float) -> FloatArray:
    """
    Quantile function of the GED(nu).
    Uses the relationship to the gamma distribution.
    """
    from scipy.stats import gamma as gamma_dist

    p = np.asarray(p, dtype=np.float64)
    lam = _ged_lambda(nu)

    # |x/lam|^nu ~ Gamma(1/nu, 2)
    # For x > 0: p(X < x) = 0.5 + 0.5 * p(Gamma(1/nu) < (x/lam)^nu)
    # Invert numerically via scipy.special
    from scipy.special import gammaincinv

    z = 2 * np.abs(p - 0.5)
    y = gammaincinv(1 / nu, z) ** (1 / nu) * lam
    return np.where(p >= 0.5, y, -y)

ged_score 

ged_score(x: FloatArray, nu: float) -> FloatArray

Analytic score d_log_f / d_nu for GED.

Used in MLE to estimate nu. The MATLAB toolbox uses finite differences.

Source code in src/mfe/distributions/ged.py 
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def ged_score(x: FloatArray, nu: float) -> FloatArray:
    """
    Analytic score d_log_f / d_nu for GED.

    Used in MLE to estimate nu. The MATLAB toolbox uses finite differences.
    """
    from scipy.special import digamma

    x = np.asarray(x, dtype=np.float64)
    lam = _ged_lambda(nu)

    # d/dnu [log c + log kernel]
    dlam_dnu = lam * (-2 / nu ** 2 * np.log(2) - 1 / nu ** 2 * (digamma(1 / nu) - digamma(3 / nu)))
    # ... (full derivation omitted for brevity, numerical fallback used in practice)
    # Approximate via finite differences for now
    h = 1e-5
    logpdf_plus = ged_logpdf(x, nu + h)
    logpdf_minus = ged_logpdf(x, nu - h)
    return (logpdf_plus - logpdf_minus) / (2 * h)

mvnorm_loglik 

mvnorm_loglik(data: FloatArray, sigma_t: FloatArray) -> float

Gaussian multivariate log-likelihood for a time-varying covariance sequence.

L = sum_t [ -K/2 * log(2pi) - 0.5 * log|Sigma_t| - 0.5 * x_t' Sigma_t^{-1} x_t ]

Parameters:

Name Type Description Default
data FloatArray
required
sigma_t (T, K, K) or (K, K) conditional covariance matrices 
  If (K, K), treated as a constant covariance
required

Returns:

Type Description
float — total log-likelihood
Source code in src/mfe/distributions/mvnorm.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def mvnorm_loglik(
    data: FloatArray,
    sigma_t: FloatArray,
) -> float:
    """
    Gaussian multivariate log-likelihood for a time-varying covariance sequence.

    L = sum_t [ -K/2 * log(2pi) - 0.5 * log|Sigma_t| - 0.5 * x_t' Sigma_t^{-1} x_t ]

    Parameters
    ----------
    data    : (T, K) return matrix (demeaned)
    sigma_t : (T, K, K) or (K, K) conditional covariance matrices
              If (K, K), treated as a constant covariance

    Returns
    -------
    float — total log-likelihood
    """
    data = np.asarray(data, dtype=np.float64)
    sigma_t = np.asarray(sigma_t, dtype=np.float64)
    T, K = data.shape

    if sigma_t.ndim == 2:
        # Constant covariance
        sigma_t = np.broadcast_to(sigma_t[None], (T, K, K))

    ll = 0.0
    const = -0.5 * K * _LOG2PI

    for t in range(T):
        S = sigma_t[t]
        sign, logdet = np.linalg.slogdet(S)
        if sign <= 0:
            return -np.inf
        try:
            quad = float(data[t] @ np.linalg.solve(S, data[t]))
        except np.linalg.LinAlgError:
            return -np.inf
        ll += const - 0.5 * (logdet + quad)

    return ll

mvnorm_loglik_t 

mvnorm_loglik_t(data_t: FloatArray, sigma: FloatArray) -> float

Single-observation Gaussian log-likelihood.

Parameters:

Name Type Description Default
data_t (K,) single observation
required
sigma FloatArray
required

Returns:

Type Description
float — log-likelihood of this observation
Source code in src/mfe/distributions/mvnorm.py 
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def mvnorm_loglik_t(
    data_t: FloatArray,
    sigma: FloatArray,
) -> float:
    """
    Single-observation Gaussian log-likelihood.

    Parameters
    ----------
    data_t : (K,) single observation
    sigma  : (K, K) covariance matrix

    Returns
    -------
    float — log-likelihood of this observation
    """
    K = len(data_t)
    sign, logdet = np.linalg.slogdet(sigma)
    if sign <= 0:
        return -np.inf
    try:
        quad = float(data_t @ np.linalg.solve(sigma, data_t))
    except np.linalg.LinAlgError:
        return -np.inf
    return -0.5 * (K * _LOG2PI + logdet + quad)

mahalanobis 

mahalanobis(data: FloatArray, sigma_t: FloatArray) -> FloatArray

Mahalanobis distances from the conditional mean.

d_t = sqrt( x_t' Sigma_t^{-1} x_t )

Parameters:

Name Type Description Default
data FloatArray
required
sigma_t (T, K, K) or (K, K)
required

Returns:

Type Description
(T,) array of Mahalanobis distances
Source code in src/mfe/distributions/mvnorm.py 
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def mahalanobis(
    data: FloatArray,
    sigma_t: FloatArray,
) -> FloatArray:
    """
    Mahalanobis distances from the conditional mean.

    d_t = sqrt( x_t' Sigma_t^{-1} x_t )

    Parameters
    ----------
    data    : (T, K)
    sigma_t : (T, K, K) or (K, K)

    Returns
    -------
    (T,) array of Mahalanobis distances
    """
    data = np.asarray(data, dtype=np.float64)
    sigma_t = np.asarray(sigma_t, dtype=np.float64)
    T, K = data.shape

    if sigma_t.ndim == 2:
        sigma_t = np.broadcast_to(sigma_t[None], (T, K, K))

    d = np.empty(T, dtype=np.float64)
    for t in range(T):
        try:
            q = float(data[t] @ np.linalg.solve(sigma_t[t], data[t]))
        except np.linalg.LinAlgError:
            q = np.nan
        d[t] = np.sqrt(max(q, 0.0))

    return d

standardize_mvn 

standardize_mvn(data: FloatArray, sigma_t: FloatArray) -> FloatArray

Extract standardized residuals: z_t = L_t^{-1} x_t where L_t L_t' = Sigma_t.

Useful for diagnostic checking: if the model is correct, z_t should be approximately i.i.d. N(0, I_K).

Parameters:

Name Type Description Default
data FloatArray
required
sigma_t (T, K, K) conditional covariance matrices
required

Returns:

Type Description
(T, K) standardized residuals
Source code in src/mfe/distributions/mvnorm.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def standardize_mvn(
    data: FloatArray,
    sigma_t: FloatArray,
) -> FloatArray:
    """
    Extract standardized residuals: z_t = L_t^{-1} x_t where L_t L_t' = Sigma_t.

    Useful for diagnostic checking: if the model is correct, z_t should be
    approximately i.i.d. N(0, I_K).

    Parameters
    ----------
    data    : (T, K) return matrix
    sigma_t : (T, K, K) conditional covariance matrices

    Returns
    -------
    (T, K) standardized residuals
    """
    data = np.asarray(data, dtype=np.float64)
    sigma_t = np.asarray(sigma_t, dtype=np.float64)
    T, K = data.shape

    if sigma_t.ndim == 2:
        sigma_t = np.broadcast_to(sigma_t[None], (T, K, K))

    z = np.empty((T, K), dtype=np.float64)
    for t in range(T):
        try:
            L = np.linalg.cholesky(sigma_t[t])
            z[t] = np.linalg.solve(L, data[t])
        except np.linalg.LinAlgError:
            z[t] = np.nan
    return z

ged

Generalized Error Distribution (GED / Power-Exponential).

Nelson, D.B. (1991): "Conditional Heteroskedasticity in Asset Returns: A New Approach", Econometrica.

PDF: f(x; nu) = nu / (2 * lambda * Gamma(1/nu)) * exp(-0.5 * |x/lambda|^nu) where lambda = (2^{-2/nu} * Gamma(1/nu) / Gamma(3/nu))^{1/2}

Special cases: nu = 1: Laplace (double exponential) nu = 2: Normal nu -> inf: Uniform

ged_logpdf 

ged_logpdf(x: FloatArray, nu: float) -> FloatArray

Log-PDF of the GED with zero mean and unit variance.

Parameters:

Name Type Description Default
x FloatArray
required
nu shape parameter (nu > 0); nu=2 is Normal
required
Source code in src/mfe/distributions/ged.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
def ged_logpdf(x: FloatArray, nu: float) -> FloatArray:
    """
    Log-PDF of the GED with zero mean and unit variance.

    Parameters
    ----------
    x  : (T,) standardized residuals
    nu : shape parameter (nu > 0); nu=2 is Normal
    """
    x = np.asarray(x, dtype=np.float64)
    lam = _ged_lambda(nu)

    log_c = np.log(nu) - np.log(2 * lam) - np.log(gamma(1 / nu))
    log_kernel = -0.5 * np.abs(x / lam) ** nu

    return log_c + log_kernel

ged_ppf 

ged_ppf(p: FloatArray, nu: float) -> FloatArray

Quantile function of the GED(nu). Uses the relationship to the gamma distribution.

Source code in src/mfe/distributions/ged.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
def ged_ppf(p: FloatArray, nu: float) -> FloatArray:
    """
    Quantile function of the GED(nu).
    Uses the relationship to the gamma distribution.
    """
    from scipy.stats import gamma as gamma_dist

    p = np.asarray(p, dtype=np.float64)
    lam = _ged_lambda(nu)

    # |x/lam|^nu ~ Gamma(1/nu, 2)
    # For x > 0: p(X < x) = 0.5 + 0.5 * p(Gamma(1/nu) < (x/lam)^nu)
    # Invert numerically via scipy.special
    from scipy.special import gammaincinv

    z = 2 * np.abs(p - 0.5)
    y = gammaincinv(1 / nu, z) ** (1 / nu) * lam
    return np.where(p >= 0.5, y, -y)

ged_score 

ged_score(x: FloatArray, nu: float) -> FloatArray

Analytic score d_log_f / d_nu for GED.

Used in MLE to estimate nu. The MATLAB toolbox uses finite differences.

Source code in src/mfe/distributions/ged.py 
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
def ged_score(x: FloatArray, nu: float) -> FloatArray:
    """
    Analytic score d_log_f / d_nu for GED.

    Used in MLE to estimate nu. The MATLAB toolbox uses finite differences.
    """
    from scipy.special import digamma

    x = np.asarray(x, dtype=np.float64)
    lam = _ged_lambda(nu)

    # d/dnu [log c + log kernel]
    dlam_dnu = lam * (-2 / nu ** 2 * np.log(2) - 1 / nu ** 2 * (digamma(1 / nu) - digamma(3 / nu)))
    # ... (full derivation omitted for brevity, numerical fallback used in practice)
    # Approximate via finite differences for now
    h = 1e-5
    logpdf_plus = ged_logpdf(x, nu + h)
    logpdf_minus = ged_logpdf(x, nu - h)
    return (logpdf_plus - logpdf_minus) / (2 * h)

mvnorm

Multivariate normal log-likelihood and related utilities.

MFE MATLAB mvnormloglik.m equivalent. Used internally in multivariate GARCH estimation but exposed here as a clean public function because it comes up constantly in empirical work.

Also provides: - mvnorm_loglik — exact Gaussian log-likelihood for a given Sigma_t sequence - mvnorm_qmle — QMLE with a fixed Sigma_t from any multivariate model - standardize_mvn — extract Mahalanobis residuals from (T, K, K) covariances

mvnorm_loglik 

mvnorm_loglik(data: FloatArray, sigma_t: FloatArray) -> float

Gaussian multivariate log-likelihood for a time-varying covariance sequence.

L = sum_t [ -K/2 * log(2pi) - 0.5 * log|Sigma_t| - 0.5 * x_t' Sigma_t^{-1} x_t ]

Parameters:

Name Type Description Default
data FloatArray
required
sigma_t (T, K, K) or (K, K) conditional covariance matrices 
  If (K, K), treated as a constant covariance
required

Returns:

Type Description
float — total log-likelihood
Source code in src/mfe/distributions/mvnorm.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def mvnorm_loglik(
    data: FloatArray,
    sigma_t: FloatArray,
) -> float:
    """
    Gaussian multivariate log-likelihood for a time-varying covariance sequence.

    L = sum_t [ -K/2 * log(2pi) - 0.5 * log|Sigma_t| - 0.5 * x_t' Sigma_t^{-1} x_t ]

    Parameters
    ----------
    data    : (T, K) return matrix (demeaned)
    sigma_t : (T, K, K) or (K, K) conditional covariance matrices
              If (K, K), treated as a constant covariance

    Returns
    -------
    float — total log-likelihood
    """
    data = np.asarray(data, dtype=np.float64)
    sigma_t = np.asarray(sigma_t, dtype=np.float64)
    T, K = data.shape

    if sigma_t.ndim == 2:
        # Constant covariance
        sigma_t = np.broadcast_to(sigma_t[None], (T, K, K))

    ll = 0.0
    const = -0.5 * K * _LOG2PI

    for t in range(T):
        S = sigma_t[t]
        sign, logdet = np.linalg.slogdet(S)
        if sign <= 0:
            return -np.inf
        try:
            quad = float(data[t] @ np.linalg.solve(S, data[t]))
        except np.linalg.LinAlgError:
            return -np.inf
        ll += const - 0.5 * (logdet + quad)

    return ll

mvnorm_loglik_t 

mvnorm_loglik_t(data_t: FloatArray, sigma: FloatArray) -> float

Single-observation Gaussian log-likelihood.

Parameters:

Name Type Description Default
data_t (K,) single observation
required
sigma FloatArray
required

Returns:

Type Description
float — log-likelihood of this observation
Source code in src/mfe/distributions/mvnorm.py 
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def mvnorm_loglik_t(
    data_t: FloatArray,
    sigma: FloatArray,
) -> float:
    """
    Single-observation Gaussian log-likelihood.

    Parameters
    ----------
    data_t : (K,) single observation
    sigma  : (K, K) covariance matrix

    Returns
    -------
    float — log-likelihood of this observation
    """
    K = len(data_t)
    sign, logdet = np.linalg.slogdet(sigma)
    if sign <= 0:
        return -np.inf
    try:
        quad = float(data_t @ np.linalg.solve(sigma, data_t))
    except np.linalg.LinAlgError:
        return -np.inf
    return -0.5 * (K * _LOG2PI + logdet + quad)

mahalanobis 

mahalanobis(data: FloatArray, sigma_t: FloatArray) -> FloatArray

Mahalanobis distances from the conditional mean.

d_t = sqrt( x_t' Sigma_t^{-1} x_t )

Parameters:

Name Type Description Default
data FloatArray
required
sigma_t (T, K, K) or (K, K)
required

Returns:

Type Description
(T,) array of Mahalanobis distances
Source code in src/mfe/distributions/mvnorm.py 
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def mahalanobis(
    data: FloatArray,
    sigma_t: FloatArray,
) -> FloatArray:
    """
    Mahalanobis distances from the conditional mean.

    d_t = sqrt( x_t' Sigma_t^{-1} x_t )

    Parameters
    ----------
    data    : (T, K)
    sigma_t : (T, K, K) or (K, K)

    Returns
    -------
    (T,) array of Mahalanobis distances
    """
    data = np.asarray(data, dtype=np.float64)
    sigma_t = np.asarray(sigma_t, dtype=np.float64)
    T, K = data.shape

    if sigma_t.ndim == 2:
        sigma_t = np.broadcast_to(sigma_t[None], (T, K, K))

    d = np.empty(T, dtype=np.float64)
    for t in range(T):
        try:
            q = float(data[t] @ np.linalg.solve(sigma_t[t], data[t]))
        except np.linalg.LinAlgError:
            q = np.nan
        d[t] = np.sqrt(max(q, 0.0))

    return d

standardize_mvn 

standardize_mvn(data: FloatArray, sigma_t: FloatArray) -> FloatArray

Extract standardized residuals: z_t = L_t^{-1} x_t where L_t L_t' = Sigma_t.

Useful for diagnostic checking: if the model is correct, z_t should be approximately i.i.d. N(0, I_K).

Parameters:

Name Type Description Default
data FloatArray
required
sigma_t (T, K, K) conditional covariance matrices
required

Returns:

Type Description
(T, K) standardized residuals
Source code in src/mfe/distributions/mvnorm.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def standardize_mvn(
    data: FloatArray,
    sigma_t: FloatArray,
) -> FloatArray:
    """
    Extract standardized residuals: z_t = L_t^{-1} x_t where L_t L_t' = Sigma_t.

    Useful for diagnostic checking: if the model is correct, z_t should be
    approximately i.i.d. N(0, I_K).

    Parameters
    ----------
    data    : (T, K) return matrix
    sigma_t : (T, K, K) conditional covariance matrices

    Returns
    -------
    (T, K) standardized residuals
    """
    data = np.asarray(data, dtype=np.float64)
    sigma_t = np.asarray(sigma_t, dtype=np.float64)
    T, K = data.shape

    if sigma_t.ndim == 2:
        sigma_t = np.broadcast_to(sigma_t[None], (T, K, K))

    z = np.empty((T, K), dtype=np.float64)
    for t in range(T):
        try:
            L = np.linalg.cholesky(sigma_t[t])
            z[t] = np.linalg.solve(L, data[t])
        except np.linalg.LinAlgError:
            z[t] = np.nan
    return z

skewt

Hansen's Skewed Student-t distribution (Hansen 1994).

Hansen, B.E. (1994): "Autoregressive Conditional Density Estimation", International Economic Review, 35(3), 705-730.

Parameters: nu : degrees of freedom (nu > 2) lam : skewness parameter (-1 < lam < 1)

PDF: f(x; nu, lam) = bc * (1 + 1/(nu-2) * ((bx+a)/(1+/-lam))^2)^{-(nu+1)/2}

with sign depending on whether x < -a/b or x >= -a/b.

The MATLAB mfe-toolbox computes gradients numerically. We provide analytic score functions for faster GARCH estimation.

skewt_logpdf 

skewt_logpdf(x: FloatArray, nu: float, lam: float) -> FloatArray

Log-PDF of Hansen's Skew-t distribution.

Parameters:

Name Type Description Default
x FloatArray
required
nu float
required
lam skewness(-1, 1)
required
Source code in src/mfe/distributions/skewt.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def skewt_logpdf(
    x: FloatArray,
    nu: float,
    lam: float,
) -> FloatArray:
    """
    Log-PDF of Hansen's Skew-t distribution.

    Parameters
    ----------
    x   : (T,) standardized residuals
    nu  : degrees of freedom (>2)
    lam : skewness (-1, 1)
    """
    x = np.asarray(x, dtype=np.float64)
    a, b, c = _hansen_constants(nu, lam)

    mask = x < -a / b
    z = np.where(mask,
                 (b * x + a) / (1 - lam),
                 (b * x + a) / (1 + lam))

    log_kernel = -(nu + 1) / 2 * np.log(1 + z ** 2 / (nu - 2))
    log_pdf = np.log(b) + np.log(c) + log_kernel

    return log_pdf

skewt_score 

skewt_score(x: FloatArray, nu: float, lam: float) -> tuple[FloatArray, FloatArray]

Analytic score of log-skewt PDF with respect to (nu, lam).

Returns (d_log_f / d_nu, d_log_f / d_lam), each (T,).

These are used in the outer MLE loop over distribution parameters, avoiding the finite-difference approach from the MATLAB source.

Source code in src/mfe/distributions/skewt.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def skewt_score(
    x: FloatArray,
    nu: float,
    lam: float,
) -> tuple[FloatArray, FloatArray]:
    """
    Analytic score of log-skewt PDF with respect to (nu, lam).

    Returns (d_log_f / d_nu, d_log_f / d_lam), each (T,).

    These are used in the outer MLE loop over distribution parameters,
    avoiding the finite-difference approach from the MATLAB source.
    """
    from scipy.special import digamma

    x = np.asarray(x, dtype=np.float64)
    a, b, c = _hansen_constants(nu, lam)

    mask = x < -a / b
    sign = np.where(mask, 1 - lam, 1 + lam)
    z = (b * x + a) / sign

    u = z ** 2 / (nu - 2)
    denom = 1 + u

    # d/d_nu: from log kernel
    d_kernel_dnu = (
        -0.5 * np.log(denom)
        + (nu + 1) / 2 * z ** 2 / ((nu - 2) ** 2 * denom)
    )
    d_c_dnu = 0.5 * (digamma((nu + 1) / 2) - digamma(nu / 2) - 1 / (nu - 2))
    score_nu = d_c_dnu + d_kernel_dnu

    # d/d_lam: via chain rule through z(lam) and a(lam), b(lam)
    # This is the expensive part — see Hansen (1994) Appendix
    da_dlam = 4 * c * (nu - 2) / (nu - 1)  # simplified: ignoring dc/dlam for now
    db_dlam = (6 * lam - 2 * a * da_dlam) / (2 * b)

    dz_dlam_pos = (db_dlam * x + da_dlam) * (1 + lam) - (b * x + a)
    dz_dlam_pos /= (1 + lam) ** 2
    dz_dlam_neg = (db_dlam * x + da_dlam) * (1 - lam) + (b * x + a)
    dz_dlam_neg /= (1 - lam) ** 2

    dz_dlam = np.where(mask, dz_dlam_neg, dz_dlam_pos)
    d_kernel_dlam = -(nu + 1) * z * dz_dlam / ((nu - 2) * denom)
    score_lam = db_dlam / b + d_kernel_dlam

    return score_nu, score_lam

skewt_ppf 

skewt_ppf(p: FloatArray, nu: float, lam: float) -> FloatArray

Quantile function (inverse CDF) of Hansen's Skew-t. Used for VaR/ES computation.

Source code in src/mfe/distributions/skewt.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def skewt_ppf(
    p: FloatArray,
    nu: float,
    lam: float,
) -> FloatArray:
    """
    Quantile function (inverse CDF) of Hansen's Skew-t.
    Used for VaR/ES computation.
    """
    p = np.asarray(p, dtype=np.float64)
    a, b, c = _hansen_constants(nu, lam)

    p1 = (1 - lam) / 2  # mass in the left tail

    result = np.where(
        p < p1,
        (student_t.ppf(p / (1 - lam), df=nu) * (1 - lam) - a) / b,
        (student_t.ppf((p - (1 - lam) / 2) / (1 + lam) + 0.5, df=nu) * (1 + lam) - a) / b,
    )
    return result