Skip to content

mfe.utils

mfe.utils

lag_matrix 

lag_matrix(x: FloatArray, lags: int | list[int], trim: bool = True) -> FloatArray

Construct a lag matrix from a 1-D or 2-D array.

Parameters:

Name Type Description Default
x array of shape (T,) or (T, K)
required
lags int or list[int]

If int, lags 1..lags are included. If list, exactly those lags are included.

 required
trim bool

If True (default), drop the leading NaN rows.

 True

Returns:

Type Description
array of shape (T - max_lag, len(lags) * K) if trim else (T, ...)
Source code in src/mfe/utils/lags.py 
15
16
17
18
19
20
21
22
23
24
25
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
def lag_matrix(x: FloatArray, lags: int | list[int], trim: bool = True) -> FloatArray:
    """
    Construct a lag matrix from a 1-D or 2-D array.

    Parameters
    ----------
    x : array of shape (T,) or (T, K)
    lags : int or list[int]
        If int, lags 1..lags are included.
        If list, exactly those lags are included.
    trim : bool
        If True (default), drop the leading NaN rows.

    Returns
    -------
    array of shape (T - max_lag, len(lags) * K) if trim else (T, ...)
    """
    x = np.asarray(x, dtype=np.float64)
    if x.ndim == 1:
        x = x[:, None]
    T, K = x.shape

    if isinstance(lags, int):
        lag_list = list(range(1, lags + 1))
    else:
        lag_list = list(lags)

    max_lag = max(lag_list)
    out = np.empty((T, len(lag_list) * K), dtype=np.float64)
    out[:] = np.nan

    for col, lag in enumerate(lag_list):
        out[lag:, col * K : (col + 1) * K] = x[: T - lag]

    if trim:
        out = out[max_lag:]

    return out

har_lag_matrix 

har_lag_matrix(rv: FloatArray, horizons: tuple[int, int, int] = (1, 5, 22), trim: bool = True) -> FloatArray

Build the HAR regressor matrix [RV_d, RV_w, RV_m] from daily RV.

Uses rolling averages, not raw lags, matching the Corsi (2009) definition: RV_{t|t-h} = (1/h) * sum_{k=0}^{h-1} RV_{t-k}

Parameters:

Name Type Description Default
rv (T,) array of daily realized variances
required
horizons tuple of 3 ints (daily, weekly, monthly)
(1, 5, 22)
trim bool — drop leading NaNs
True

Returns:

Type Description
(T - max_h, 3) array: columns are [RV_d_lag1, RV_w, RV_m]
Source code in src/mfe/utils/lags.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
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
95
96
def har_lag_matrix(
    rv: FloatArray,
    horizons: tuple[int, int, int] = (1, 5, 22),
    trim: bool = True,
) -> FloatArray:
    """
    Build the HAR regressor matrix [RV_d, RV_w, RV_m] from daily RV.

    Uses rolling averages, not raw lags, matching the Corsi (2009) definition:
        RV_{t|t-h} = (1/h) * sum_{k=0}^{h-1} RV_{t-k}

    Parameters
    ----------
    rv : (T,) array of daily realized variances
    horizons : tuple of 3 ints (daily, weekly, monthly), default (1, 5, 22)
    trim : bool — drop leading NaNs

    Returns
    -------
    (T - max_h, 3) array: columns are [RV_d_lag1, RV_w, RV_m]
    """
    rv = np.asarray(rv, dtype=np.float64)
    T = len(rv)
    max_h = max(horizons)

    cols = []
    for h in horizons:
        # rolling mean over h observations, shifted by 1 (yesterday's average)
        kernel = np.ones(h) / h
        rolled = np.convolve(rv, kernel, mode="full")[:T]
        # shift right by 1: today's regressor is yesterday's rolling mean
        col = np.empty(T, dtype=np.float64)
        col[:] = np.nan
        col[1:] = rolled[:-1]
        cols.append(col)

    out = np.column_stack(cols)

    if trim:
        out = out[max_h:]

    return out

sandwich 

sandwich(scores: FloatArray, hessian: FloatArray) -> FloatArray

Sandwich (robust) covariance: H^{-1} B H^{-1} where B = scores.T @ scores / T and H is the (negative) Hessian / T.

Parameters:

Name Type Description Default
scores FloatArray
required
hessian (P, P) negative Hessian evaluated at MLE
required

Returns:

Type Description
(P, P) robust covariance matrix
Source code in src/mfe/utils/vcv.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def sandwich(
    scores: FloatArray,
    hessian: FloatArray,
) -> FloatArray:
    """
    Sandwich (robust) covariance: H^{-1} B H^{-1}
    where B = scores.T @ scores / T and H is the (negative) Hessian / T.

    Parameters
    ----------
    scores  : (T, P) score matrix
    hessian : (P, P) negative Hessian evaluated at MLE

    Returns
    -------
    (P, P) robust covariance matrix
    """
    T = scores.shape[0]
    B = scores.T @ scores / T
    H_inv = np.linalg.inv(hessian / T)
    return H_inv @ B @ H_inv / T

newey_west 

newey_west(scores: FloatArray, bandwidth: int | None = None, hessian: FloatArray | None = None) -> FloatArray

Newey-West (HAC) covariance.

Parameters:

Name Type Description Default
scores FloatArray
required
bandwidth number of lags; if None uses Andrews (1991) automatic selector
None
hessian FloatArray | None
None

Returns:

Type Description
(P, P) matrix: B_hat (HAC meat) if hessian is None, else full sandwich
Source code in src/mfe/utils/vcv.py 
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
68
69
70
71
72
73
def newey_west(
    scores: FloatArray,
    bandwidth: int | None = None,
    hessian: FloatArray | None = None,
) -> FloatArray:
    """
    Newey-West (HAC) covariance.

    Parameters
    ----------
    scores    : (T, P) score matrix
    bandwidth : number of lags; if None uses Andrews (1991) automatic selector
    hessian   : (P, P) — if provided, returns sandwich; otherwise returns B_hat only

    Returns
    -------
    (P, P) matrix: B_hat (HAC meat) if hessian is None, else full sandwich
    """
    T, P = scores.shape

    if bandwidth is None:
        bandwidth = int(np.floor(4 * (T / 100) ** (2 / 9)))

    # Newey-West weights: 1 - h/(bandwidth+1)
    B = scores.T @ scores / T
    for lag in range(1, bandwidth + 1):
        w = 1.0 - lag / (bandwidth + 1)
        gamma = scores[lag:].T @ scores[:T - lag] / T
        B += w * (gamma + gamma.T)

    if hessian is None:
        return B

    H_inv = np.linalg.inv(hessian / T)
    return H_inv @ B @ H_inv / T

lags

Vectorized lag-matrix utilities.

All functions operate on (T,) or (T, K) arrays and return views or stride-tricks arrays where possible — no unnecessary copies.

lag_matrix 

lag_matrix(x: FloatArray, lags: int | list[int], trim: bool = True) -> FloatArray

Construct a lag matrix from a 1-D or 2-D array.

Parameters:

Name Type Description Default
x array of shape (T,) or (T, K)
required
lags int or list[int]

If int, lags 1..lags are included. If list, exactly those lags are included.

 required
trim bool

If True (default), drop the leading NaN rows.

 True

Returns:

Type Description
array of shape (T - max_lag, len(lags) * K) if trim else (T, ...)
Source code in src/mfe/utils/lags.py 
15
16
17
18
19
20
21
22
23
24
25
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
def lag_matrix(x: FloatArray, lags: int | list[int], trim: bool = True) -> FloatArray:
    """
    Construct a lag matrix from a 1-D or 2-D array.

    Parameters
    ----------
    x : array of shape (T,) or (T, K)
    lags : int or list[int]
        If int, lags 1..lags are included.
        If list, exactly those lags are included.
    trim : bool
        If True (default), drop the leading NaN rows.

    Returns
    -------
    array of shape (T - max_lag, len(lags) * K) if trim else (T, ...)
    """
    x = np.asarray(x, dtype=np.float64)
    if x.ndim == 1:
        x = x[:, None]
    T, K = x.shape

    if isinstance(lags, int):
        lag_list = list(range(1, lags + 1))
    else:
        lag_list = list(lags)

    max_lag = max(lag_list)
    out = np.empty((T, len(lag_list) * K), dtype=np.float64)
    out[:] = np.nan

    for col, lag in enumerate(lag_list):
        out[lag:, col * K : (col + 1) * K] = x[: T - lag]

    if trim:
        out = out[max_lag:]

    return out

har_lag_matrix 

har_lag_matrix(rv: FloatArray, horizons: tuple[int, int, int] = (1, 5, 22), trim: bool = True) -> FloatArray

Build the HAR regressor matrix [RV_d, RV_w, RV_m] from daily RV.

Uses rolling averages, not raw lags, matching the Corsi (2009) definition: RV_{t|t-h} = (1/h) * sum_{k=0}^{h-1} RV_{t-k}

Parameters:

Name Type Description Default
rv (T,) array of daily realized variances
required
horizons tuple of 3 ints (daily, weekly, monthly)
(1, 5, 22)
trim bool — drop leading NaNs
True

Returns:

Type Description
(T - max_h, 3) array: columns are [RV_d_lag1, RV_w, RV_m]
Source code in src/mfe/utils/lags.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
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
95
96
def har_lag_matrix(
    rv: FloatArray,
    horizons: tuple[int, int, int] = (1, 5, 22),
    trim: bool = True,
) -> FloatArray:
    """
    Build the HAR regressor matrix [RV_d, RV_w, RV_m] from daily RV.

    Uses rolling averages, not raw lags, matching the Corsi (2009) definition:
        RV_{t|t-h} = (1/h) * sum_{k=0}^{h-1} RV_{t-k}

    Parameters
    ----------
    rv : (T,) array of daily realized variances
    horizons : tuple of 3 ints (daily, weekly, monthly), default (1, 5, 22)
    trim : bool — drop leading NaNs

    Returns
    -------
    (T - max_h, 3) array: columns are [RV_d_lag1, RV_w, RV_m]
    """
    rv = np.asarray(rv, dtype=np.float64)
    T = len(rv)
    max_h = max(horizons)

    cols = []
    for h in horizons:
        # rolling mean over h observations, shifted by 1 (yesterday's average)
        kernel = np.ones(h) / h
        rolled = np.convolve(rv, kernel, mode="full")[:T]
        # shift right by 1: today's regressor is yesterday's rolling mean
        col = np.empty(T, dtype=np.float64)
        col[:] = np.nan
        col[1:] = rolled[:-1]
        cols.append(col)

    out = np.column_stack(cols)

    if trim:
        out = out[max_h:]

    return out

typing

Shared type aliases for the mfe package.

vcv

Robust covariance matrix estimators.

  • sandwich (QMLE robust, a.k.a. Huber-White)
  • newey_west (HAC)
  • Both accept a (T, P) score matrix and optionally a (P, P) Hessian.

sandwich 

sandwich(scores: FloatArray, hessian: FloatArray) -> FloatArray

Sandwich (robust) covariance: H^{-1} B H^{-1} where B = scores.T @ scores / T and H is the (negative) Hessian / T.

Parameters:

Name Type Description Default
scores FloatArray
required
hessian (P, P) negative Hessian evaluated at MLE
required

Returns:

Type Description
(P, P) robust covariance matrix
Source code in src/mfe/utils/vcv.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def sandwich(
    scores: FloatArray,
    hessian: FloatArray,
) -> FloatArray:
    """
    Sandwich (robust) covariance: H^{-1} B H^{-1}
    where B = scores.T @ scores / T and H is the (negative) Hessian / T.

    Parameters
    ----------
    scores  : (T, P) score matrix
    hessian : (P, P) negative Hessian evaluated at MLE

    Returns
    -------
    (P, P) robust covariance matrix
    """
    T = scores.shape[0]
    B = scores.T @ scores / T
    H_inv = np.linalg.inv(hessian / T)
    return H_inv @ B @ H_inv / T

newey_west 

newey_west(scores: FloatArray, bandwidth: int | None = None, hessian: FloatArray | None = None) -> FloatArray

Newey-West (HAC) covariance.

Parameters:

Name Type Description Default
scores FloatArray
required
bandwidth number of lags; if None uses Andrews (1991) automatic selector
None
hessian FloatArray | None
None

Returns:

Type Description
(P, P) matrix: B_hat (HAC meat) if hessian is None, else full sandwich
Source code in src/mfe/utils/vcv.py 
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
68
69
70
71
72
73
def newey_west(
    scores: FloatArray,
    bandwidth: int | None = None,
    hessian: FloatArray | None = None,
) -> FloatArray:
    """
    Newey-West (HAC) covariance.

    Parameters
    ----------
    scores    : (T, P) score matrix
    bandwidth : number of lags; if None uses Andrews (1991) automatic selector
    hessian   : (P, P) — if provided, returns sandwich; otherwise returns B_hat only

    Returns
    -------
    (P, P) matrix: B_hat (HAC meat) if hessian is None, else full sandwich
    """
    T, P = scores.shape

    if bandwidth is None:
        bandwidth = int(np.floor(4 * (T / 100) ** (2 / 9)))

    # Newey-West weights: 1 - h/(bandwidth+1)
    B = scores.T @ scores / T
    for lag in range(1, bandwidth + 1):
        w = 1.0 - lag / (bandwidth + 1)
        gamma = scores[lag:].T @ scores[:T - lag] / T
        B += w * (gamma + gamma.T)

    if hessian is None:
        return B

    H_inv = np.linalg.inv(hessian / T)
    return H_inv @ B @ H_inv / T