mfe.timeseries¶

mfe.timeseries ¶

mfe.timeseries — Time series models.

vectorar VAR(P) estimation with 4 VCV options grangercause Granger causality LR/LM/Wald tests impulse_response IRF with delta-method standard errors beveridge_nelson Beveridge-Nelson trend/cycle decomposition for I(1) series

VARResult `dataclass` ¶

VARResult(params: list[FloatArray], const: FloatArray | None, errors: FloatArray, sigma: FloatArray, r_squared: FloatArray, vcv: FloatArray, param_vec: FloatArray, lags: list[int], n_obs: int, n_vars: int, log_likelihood: float)

VAR(P) estimation result.

GCResult `dataclass` ¶

GCResult(statistics: FloatArray, p_values: FloatArray, method: str, n_obs: int, n_vars: int)

Granger causality test result.

IRFResult `dataclass` ¶

IRFResult(responses: FloatArray, std_errors: FloatArray, lags: int, horizon: int, decomp: str)

Impulse response function result.

BNResult `dataclass` ¶

BNResult(trend: FloatArray, cycle: FloatArray, original: FloatArray, drift: float, ar_params: FloatArray, ar_order: int, method: str)

Beveridge-Nelson decomposition result.

vectorar ¶

vectorar(y: FloatArray, lags: int | list[int] = 1, include_const: bool = True, het: bool = True, uncorr: bool = False) -> VARResult

Estimate a VAR(P) or irregular VAR.

Parameters:

Name	Type	Default
`y`	`FloatArray`	required
`lags`	`int \| list[int]`	`1`
`include_const`	`bool`	`True`
`het`	`bool`	`True`
`uncorr`	`bool`	`False`

Returns:

Type	Description
`VARResult with params as list of (K,K) matrices, one per lag.`

Source code in src/mfe/timeseries/var.py

def vectorar(
    y: FloatArray,
    lags: int | list[int] = 1,
    include_const: bool = True,
    het: bool = True,
    uncorr: bool = False,
) -> VARResult:
    """
    Estimate a VAR(P) or irregular VAR.

    Parameters
    ----------
    y            : (T, K) data matrix
    lags         : int (regular VAR(P)) or list[int] (irregular VAR, e.g. [1,3])
    include_const: include a constant term (default True)
    het          : heteroskedasticity-robust VCV (default True)
    uncorr       : assume uncorrelated errors across equations (default False)

    Returns
    -------
    VARResult with params as list of (K,K) matrices, one per lag.
    """
    Y = np.asarray(y, dtype=np.float64)
    T, K = Y.shape

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

    Y_adj, X, max_lag = _build_regressor_matrix(Y, lag_list, include_const)
    n, n_regs = X.shape
    n_lag_params = len(lag_list) * K

    # OLS equation-by-equation (GLS = OLS for VAR)
    try:
        XTX_inv = np.linalg.inv(X.T @ X)
    except np.linalg.LinAlgError:
        raise ValueError("Regressors are singular — reduce lag order or check data.")

    B_full = XTX_inv @ X.T @ Y_adj   # (n_regs, K)
    errors = Y_adj - X @ B_full       # (n, K)
    Sigma = errors.T @ errors / n

    # Unpack constant and lag matrices
    offset = 0
    const_vec = None
    if include_const:
        const_vec = B_full[0]       # (K,)
        offset = 1

    params = []
    for lag in lag_list:
        block = B_full[offset: offset + K]   # (K, K) — rows = regressors for this lag
        params.append(block.T)               # (K, K) as Phi_p (rows = equations)
        offset += K

    # Build param_vec in the MFE ordering:
    # for each equation k: [const_k (opt), phi_{k,1,1}..phi_{k,1,K}, phi_{k,2,1}..., ...]
    param_parts = []
    for k in range(K):
        row = []
        if include_const:
            row.append([const_vec[k]])
        for p_idx in range(len(lag_list)):
            row.append(params[p_idx][k])
        param_parts.append(np.concatenate(row))
    param_vec = np.concatenate(param_parts)

    # VCV
    if het and not uncorr:
        vcv = _vcv_het_corr(X, errors, K, n_regs)
    elif het and uncorr:
        vcv = _vcv_het_uncorr(X, errors, K, n_regs)
    elif not het and uncorr:
        vcv = _vcv_hom_uncorr(X, errors, K, n_regs)
    else:
        vcv = _vcv_hom_corr(X, errors, K, n_regs)

    # R^2 per equation
    r2 = np.array([
        1.0 - float(np.sum(errors[:, k] ** 2)) / float(np.sum((Y_adj[:, k] - Y_adj[:, k].mean()) ** 2))
        for k in range(K)
    ])

    # Log-likelihood (Gaussian)
    sign, logdet = np.linalg.slogdet(Sigma)
    ll = -0.5 * n * (K * np.log(2 * np.pi) + logdet + K) if sign > 0 else -1e10

    return VARResult(
        params=params,
        const=const_vec,
        errors=errors,
        sigma=Sigma,
        r_squared=r2,
        vcv=vcv,
        param_vec=param_vec,
        lags=lag_list,
        n_obs=n,
        n_vars=K,
        log_likelihood=float(ll),
    )

grangercause ¶

grangercause(y: FloatArray, lags: int | list[int] = 1, include_const: bool = True, het: bool = True, uncorr: bool = False, method: Literal['lr', 'lm', 'wald'] = 'lr') -> GCResult

Granger causality testing in a VAR.

stat[i,j] tests H0: lags of y_j do not Granger-cause y_i.

Parameters:

Name	Type	Description	Default
`y`	`same as vectorar`		required
`lags`	`same as vectorar`		required
`include_const`	`same as vectorar`		required
`het`	`same as vectorar`		required
`uncorr`	`same as vectorar`		required
`method`	`'lr' \| 'lm' \| 'wald'`	"lr" — likelihood ratio (Chi2, robust if het=True) "lm" — score/LM test "wald" — Wald test using VCV from vectorar	`'lr'`

Returns:

Type	Description
`GCResult with (K,K) matrices of statistics and p-values.`

Source code in src/mfe/timeseries/var.py

def grangercause(
    y: FloatArray,
    lags: int | list[int] = 1,
    include_const: bool = True,
    het: bool = True,
    uncorr: bool = False,
    method: Literal["lr", "lm", "wald"] = "lr",
) -> GCResult:
    """
    Granger causality testing in a VAR.

    stat[i,j] tests H0: lags of y_j do not Granger-cause y_i.

    Parameters
    ----------
    y, lags, include_const, het, uncorr : same as vectorar
    method : "lr" | "lm" | "wald"
        "lr"   — likelihood ratio (Chi2, robust if het=True)
        "lm"   — score/LM test
        "wald" — Wald test using VCV from vectorar

    Returns
    -------
    GCResult with (K,K) matrices of statistics and p-values.
    """
    Y = np.asarray(y, dtype=np.float64)
    T, K = Y.shape

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

    P = len(lag_list)
    df = P  # number of restrictions per (i,j) pair

    res_unr = vectorar(Y, lags=lag_list, include_const=include_const, het=het, uncorr=uncorr)

    stat_mat = np.full((K, K), np.nan, dtype=np.float64)
    pval_mat = np.full((K, K), np.nan, dtype=np.float64)

    Y_adj, X_full, max_lag = _build_regressor_matrix(Y, lag_list, include_const)
    n = Y_adj.shape[0]
    n_regs = X_full.shape[1]
    n_const = 1 if include_const else 0

    for i in range(K):       # caused variable (equation)
        for j in range(K):   # causing variable (excluded)
            if i == j:
                stat_mat[i, j] = np.nan
                pval_mat[i, j] = np.nan
                continue

            # Build restricted X: drop columns for lags of y_j
            # Column layout: [const?] [lag1: K cols] [lag2: K cols] ...
            drop_cols = []
            for p in range(P):
                col_start = n_const + p * K
                drop_cols.append(col_start + j)   # j-th variable in lag p block
            keep_cols = [c for c in range(n_regs) if c not in drop_cols]
            X_r = X_full[:, keep_cols]

            y_i = Y_adj[:, i]

            if method == "wald":
                # Wald: R beta = 0 using VCV from unrestricted
                # Extract VCV block for equation i
                eqs_per_col = n_regs
                vcv_i = res_unr.vcv[i*eqs_per_col:(i+1)*eqs_per_col,
                                     i*eqs_per_col:(i+1)*eqs_per_col]
                # Restriction matrix R: picks rows corresponding to j's lags
                R = np.zeros((df, n_regs), dtype=np.float64)
                for p_idx, col in enumerate(drop_cols):
                    R[p_idx, col] = 1.0
                # Param vector for equation i
                B_i = res_unr.param_vec[i*n_regs:(i+1)*n_regs]
                Rb = R @ B_i
                try:
                    W = float(n * Rb @ np.linalg.solve(R @ vcv_i @ R.T * n, Rb))
                except np.linalg.LinAlgError:
                    W = np.nan
                stat_mat[i, j] = W
                pval_mat[i, j] = float(1 - stats.chi2.cdf(W, df=df)) if np.isfinite(W) else np.nan

            else:
                # LR or LM: compare restricted and unrestricted models for equation i
                # Unrestricted
                XTX_inv_u = np.linalg.inv(X_full.T @ X_full)
                b_u = XTX_inv_u @ (X_full.T @ y_i)
                e_u = y_i - X_full @ b_u
                s2_u = float(e_u @ e_u) / n

                # Restricted
                try:
                    XTX_inv_r = np.linalg.inv(X_r.T @ X_r)
                except np.linalg.LinAlgError:
                    continue
                b_r = XTX_inv_r @ (X_r.T @ y_i)
                e_r = y_i - X_r @ b_r
                s2_r = float(e_r @ e_r) / n

                if method == "lr":
                    if het:
                        # Robust LR-class: based on scores under null, VCV under alt
                        scores = X_full[:, drop_cols] * e_u[:, None]
                        B_meat = scores.T @ scores / n
                        B_bread = X_full[:, drop_cols].T @ X_full[:, drop_cols] / n
                        try:
                            B_inv = np.linalg.inv(B_bread)
                            S_inv = np.linalg.inv(B_inv @ B_meat @ B_inv / n)
                        except np.linalg.LinAlgError:
                            continue
                        s_bar = scores.mean(axis=0)
                        LR = float(n * s_bar @ S_inv @ s_bar)
                    else:
                        LR = float(n * (np.log(s2_r) - np.log(s2_u)))
                    stat_mat[i, j] = LR
                    pval_mat[i, j] = float(1 - stats.chi2.cdf(LR, df=df))

                else:  # lm
                    # LM: regress e_r on X_full, R^2 * n
                    b_aux = np.linalg.lstsq(X_full, e_r, rcond=None)[0]
                    e_aux = e_r - X_full @ b_aux
                    ss_res_aux = float(e_aux @ e_aux)
                    ss_tot_aux = float(e_r @ e_r)
                    r2_aux = 1.0 - ss_res_aux / ss_tot_aux if ss_tot_aux > 0 else 0.0
                    LM = float(n * r2_aux)
                    stat_mat[i, j] = LM
                    pval_mat[i, j] = float(1 - stats.chi2.cdf(LM, df=df))

    return GCResult(
        statistics=stat_mat,
        p_values=pval_mat,
        method=method,
        n_obs=n,
        n_vars=K,
    )

impulse_response ¶

impulse_response(y: FloatArray, lags: int | list[int] = 1, horizon: int = 12, decomp: Literal['unit', 'cholesky', 'spectral'] = 'cholesky', include_const: bool = True, het: bool = True, uncorr: bool = False) -> IRFResult

Impulse response functions for a VAR(P) with delta-method standard errors.

Parameters:

Name	Type	Description	Default
`y`	`FloatArray`		required
`lags`	`int \| list[int]`		`1`
`horizon`	`int`		`12`
`decomp`	`Literal['unit', 'cholesky', 'spectral']`	"unit" — unit shocks (unscaled), i.e. P0 = I_K "cholesky" — Cholesky of Sigma (lower triangular), recursive identification "spectral" — symmetric square root of Sigma (spectral decomposition)	`'cholesky'`
`het`	`VCV options for standard error computation`		`True`
`uncorr`	`VCV options for standard error computation`		`True`

Returns:

Type	Description
`IRFResult`	.responses : (K, K, H+1) — responses[response_var, shock_var, h] .std_errors : (K, K, H+1) delta-method std errors

Source code in src/mfe/timeseries/var.py

def impulse_response(
    y: FloatArray,
    lags: int | list[int] = 1,
    horizon: int = 12,
    decomp: Literal["unit", "cholesky", "spectral"] = "cholesky",
    include_const: bool = True,
    het: bool = True,
    uncorr: bool = False,
) -> IRFResult:
    """
    Impulse response functions for a VAR(P) with delta-method standard errors.

    Parameters
    ----------
    y        : (T, K) data
    lags     : VAR lag order or list
    horizon  : number of periods H; returns H+1 responses (including period 0)
    decomp   : shock decomposition
        "unit"      — unit shocks (unscaled), i.e. P0 = I_K
        "cholesky"  — Cholesky of Sigma (lower triangular), recursive identification
        "spectral"  — symmetric square root of Sigma (spectral decomposition)
    het, uncorr : VCV options for standard error computation

    Returns
    -------
    IRFResult
        .responses  : (K, K, H+1) — responses[response_var, shock_var, h]
        .std_errors : (K, K, H+1) delta-method std errors
    """
    Y = np.asarray(y, dtype=np.float64)
    T, K = Y.shape

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

    res = vectorar(Y, lags=lag_list, include_const=include_const, het=het, uncorr=uncorr)
    P = len(lag_list)
    Sigma = res.sigma

    # Shock decomposition matrix P0
    if decomp == "unit":
        P0 = np.eye(K)
    elif decomp == "cholesky":
        try:
            P0 = np.linalg.cholesky(Sigma)    # lower triangular
        except np.linalg.LinAlgError:
            warnings.warn("Cholesky failed; using spectral decomposition.", RuntimeWarning)
            P0 = _spectral_sqrt(Sigma)
    elif decomp == "spectral":
        P0 = _spectral_sqrt(Sigma)
    else:
        raise ValueError(f"decomp must be 'unit', 'cholesky', or 'spectral', got '{decomp}'")

    # Companion form: convert VAR(P) to VAR(1) state A of size (P*K, P*K)
    # A = [[Phi_1, Phi_2, ..., Phi_P],
    #      [I_K,   0,    ..., 0     ],
    #      [0,     I_K,  ..., 0     ],
    #      ...                       ]
    max_lag = max(lag_list)
    A = np.zeros((max_lag * K, max_lag * K), dtype=np.float64)

    for p_idx, p in enumerate(lag_list):
        # res.params[p_idx] is (K, K) Phi_p
        A[:K, p_idx * K: (p_idx + 1) * K] = res.params[p_idx]
    for block in range(1, max_lag):
        A[block * K: (block + 1) * K, (block - 1) * K: block * K] = np.eye(K)

    # IRF via MA(inf) representation: Psi_h = J A^h J' where J = [I_K | 0]
    J = np.zeros((K, max_lag * K), dtype=np.float64)
    J[:K, :K] = np.eye(K)

    responses  = np.empty((K, K, horizon + 1), dtype=np.float64)
    std_errors = np.empty((K, K, horizon + 1), dtype=np.float64)

    Ah = np.eye(max_lag * K, dtype=np.float64)
    for h in range(horizon + 1):
        Psi_h = J @ Ah @ J.T        # (K, K)
        responses[:, :, h] = Psi_h @ P0
        Ah = Ah @ A

    # Delta-method std errors via asymptotic approximation
    # For each horizon h, Var(vec(Psi_h * P0)) via chain rule on A^h
    # We use numerical differentiation for generality
    std_errors = _irf_std_errors(res, lag_list, horizon, P0, J, A, K, max_lag)

    return IRFResult(
        responses=responses,
        std_errors=std_errors,
        lags=len(lag_list),
        horizon=horizon,
        decomp=decomp,
    )

beveridge_nelson ¶

Beveridge-Nelson Decomposition.

Beveridge, S. & Nelson, C.R. (1981): "A New Approach to Decomposition of Economic Time Series into Permanent and Transitory Components with Particular Attention to Measurement of the Business Cycle", Journal of Monetary Economics, 7(2), 151-174.

The BN decomposition splits an I(1) series y_t into: y_t = tau_t + c_t

where: tau_t = permanent (trend) component — a random walk with drift c_t = transitory (cycle) component — a zero-mean stationary process

The trend is defined as the long-run forecast: tau_t = lim_{h→∞} E[y_{t+h} - h*mu | I_t]

where mu = drift of y_t = E[Delta y_t].

The cycle is: c_t = y_t - tau_t = -sum_{j=1}^{∞} E[Delta y_{t+j} - mu | I_t]

Computation

Given a forecasting model for Delta y_t (typically an AR or ARMA), the BN decomposition can be computed exactly without truncating infinite sums.

Two approaches are implemented:

State-space (exact): Cast the ARMA model into companion form and compute the long-run forecast analytically using the matrix (I - A)^{-1}. This matches the algorithm of Morley (2002) and the MFE MATLAB implementation.
Direct AR: Fit an AR(p) to Delta y by OLS and compute the BN trend via the standard formula c_t = -sum_{k=1}^{p} pi_k * Delta y_{t-k+1} (Stock & Watson 1988; Cogley 2001). Faster and simpler.

The MFE MATLAB beveridgenelson.m uses approach 2 (AR on first differences). We implement both and default to approach 2.

References

Morley, J.C. (2002): "A State–Space Approach to Calculating the Beveridge–Nelson Decomposition", Economics Letters, 75(1), 123-127.

Newbold, P. (1990): "Precise and Efficient Computation of the Beveridge–Nelson Decomposition of Economic Time Series", Journal of Monetary Economics.

BNResult `dataclass` ¶

BNResult(trend: FloatArray, cycle: FloatArray, original: FloatArray, drift: float, ar_params: FloatArray, ar_order: int, method: str)

Beveridge-Nelson decomposition result.

beveridge_nelson ¶

beveridge_nelson(y: FloatArray, ar_order: int | None = None, method: str = 'ar', ic: str = 'aic') -> BNResult

Beveridge-Nelson decomposition of an I(1) time series.

Parameters:

Name	Type	Description	Default
`y`	`FloatArray`		required
`ar_order`	`AR order p for the model of Delta y.`	`If None, selects automatically by AIC/BIC up to min(T//4, 24).`	`None`
`method`	`str`	`"state_space" — exact via companion form (Morley 2002)`	`'ar'`
`ic`	`str`		`'aic'`

Returns:

Type	Description
`BNResult`	.trend — permanent component tau_t (same length as y) .cycle — transitory component c_t = y_t - tau_t .drift — estimated drift of Delta y

Notes

The BN trend is NOT smooth — it inherits all the innovation variance of the series. If you want a smooth trend, use HP or BK filter instead. The BN cycle is zero-mean, stationary, and reflects the business-cycle component as defined by forecasts.

Source code in src/mfe/timeseries/beveridge_nelson.py

def beveridge_nelson(
    y: FloatArray,
    ar_order: int | None = None,
    method: str = "ar",
    ic: str = "aic",
) -> BNResult:
    """
    Beveridge-Nelson decomposition of an I(1) time series.

    Parameters
    ----------
    y        : (T,) level series (must be I(1) — i.e. Delta y should be stationary)
    ar_order : AR order p for the model of Delta y.
               If None, selects automatically by AIC/BIC up to min(T//4, 24).
    method   : "ar" (default) — direct AR on first differences (Cogley 2001)
               "state_space" — exact via companion form (Morley 2002)
    ic       : "aic" | "bic" — information criterion for automatic order selection

    Returns
    -------
    BNResult
        .trend  — permanent component tau_t (same length as y)
        .cycle  — transitory component c_t = y_t - tau_t
        .drift  — estimated drift of Delta y

    Notes
    -----
    The BN trend is NOT smooth — it inherits all the innovation variance of
    the series. If you want a smooth trend, use HP or BK filter instead.
    The BN cycle is zero-mean, stationary, and reflects the business-cycle
    component as defined by forecasts.
    """
    y = np.asarray(y, dtype=np.float64)
    T = len(y)
    dy = np.diff(y)               # first differences
    mu = float(dy.mean())         # drift

    # Determine AR order
    max_p = min(T // 4, 24)
    if ar_order is None:
        ar_order = _select_ar_order(dy, max_p=max_p, ic=ic)

    p = ar_order

    if method == "ar":
        trend, cycle = _bn_ar(y, dy, mu, p)
    elif method == "state_space":
        trend, cycle = _bn_state_space(y, dy, mu, p)
    else:
        raise ValueError(f"method must be 'ar' or 'state_space', got '{method}'")

    phi = _fit_ar(dy - mu, p) if p > 0 else np.array([])

    return BNResult(
        trend=trend,
        cycle=cycle,
        original=y,
        drift=mu,
        ar_params=phi,
        ar_order=p,
        method=method,
    )

var ¶

Vector Autoregression: estimation, Granger causality, impulse response functions.

vectorar — VAR(P) estimation with 4 VCV options grangercause — Granger causality LR / LM / Wald tests impulse_response — IRF with delta-method std errors, Cholesky or spectral decomp

Design gap vs. statsmodels.tsa.VAR: statsmodels VAR: OLS only, homoskedastic VCV, no heteroskedastic sandwich statsmodels IRF: standard errors only under homoskedastic assumption No Granger causality test with robust VCV in statsmodels

This module fills those gaps exactly, matching the MFE MATLAB vectorar.m outputs.

VARResult `dataclass` ¶

VARResult(params: list[FloatArray], const: FloatArray | None, errors: FloatArray, sigma: FloatArray, r_squared: FloatArray, vcv: FloatArray, param_vec: FloatArray, lags: list[int], n_obs: int, n_vars: int, log_likelihood: float)

VAR(P) estimation result.

GCResult `dataclass` ¶

GCResult(statistics: FloatArray, p_values: FloatArray, method: str, n_obs: int, n_vars: int)

Granger causality test result.

IRFResult `dataclass` ¶

IRFResult(responses: FloatArray, std_errors: FloatArray, lags: int, horizon: int, decomp: str)

Impulse response function result.

vectorar ¶

vectorar(y: FloatArray, lags: int | list[int] = 1, include_const: bool = True, het: bool = True, uncorr: bool = False) -> VARResult

Estimate a VAR(P) or irregular VAR.

Parameters:

Name	Type	Default
`y`	`FloatArray`	required
`lags`	`int \| list[int]`	`1`
`include_const`	`bool`	`True`
`het`	`bool`	`True`
`uncorr`	`bool`	`False`

Returns:

Type	Description
`VARResult with params as list of (K,K) matrices, one per lag.`

Source code in src/mfe/timeseries/var.py

def vectorar(
    y: FloatArray,
    lags: int | list[int] = 1,
    include_const: bool = True,
    het: bool = True,
    uncorr: bool = False,
) -> VARResult:
    """
    Estimate a VAR(P) or irregular VAR.

    Parameters
    ----------
    y            : (T, K) data matrix
    lags         : int (regular VAR(P)) or list[int] (irregular VAR, e.g. [1,3])
    include_const: include a constant term (default True)
    het          : heteroskedasticity-robust VCV (default True)
    uncorr       : assume uncorrelated errors across equations (default False)

    Returns
    -------
    VARResult with params as list of (K,K) matrices, one per lag.
    """
    Y = np.asarray(y, dtype=np.float64)
    T, K = Y.shape

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

    Y_adj, X, max_lag = _build_regressor_matrix(Y, lag_list, include_const)
    n, n_regs = X.shape
    n_lag_params = len(lag_list) * K

    # OLS equation-by-equation (GLS = OLS for VAR)
    try:
        XTX_inv = np.linalg.inv(X.T @ X)
    except np.linalg.LinAlgError:
        raise ValueError("Regressors are singular — reduce lag order or check data.")

    B_full = XTX_inv @ X.T @ Y_adj   # (n_regs, K)
    errors = Y_adj - X @ B_full       # (n, K)
    Sigma = errors.T @ errors / n

    # Unpack constant and lag matrices
    offset = 0
    const_vec = None
    if include_const:
        const_vec = B_full[0]       # (K,)
        offset = 1

    params = []
    for lag in lag_list:
        block = B_full[offset: offset + K]   # (K, K) — rows = regressors for this lag
        params.append(block.T)               # (K, K) as Phi_p (rows = equations)
        offset += K

    # Build param_vec in the MFE ordering:
    # for each equation k: [const_k (opt), phi_{k,1,1}..phi_{k,1,K}, phi_{k,2,1}..., ...]
    param_parts = []
    for k in range(K):
        row = []
        if include_const:
            row.append([const_vec[k]])
        for p_idx in range(len(lag_list)):
            row.append(params[p_idx][k])
        param_parts.append(np.concatenate(row))
    param_vec = np.concatenate(param_parts)

    # VCV
    if het and not uncorr:
        vcv = _vcv_het_corr(X, errors, K, n_regs)
    elif het and uncorr:
        vcv = _vcv_het_uncorr(X, errors, K, n_regs)
    elif not het and uncorr:
        vcv = _vcv_hom_uncorr(X, errors, K, n_regs)
    else:
        vcv = _vcv_hom_corr(X, errors, K, n_regs)

    # R^2 per equation
    r2 = np.array([
        1.0 - float(np.sum(errors[:, k] ** 2)) / float(np.sum((Y_adj[:, k] - Y_adj[:, k].mean()) ** 2))
        for k in range(K)
    ])

    # Log-likelihood (Gaussian)
    sign, logdet = np.linalg.slogdet(Sigma)
    ll = -0.5 * n * (K * np.log(2 * np.pi) + logdet + K) if sign > 0 else -1e10

    return VARResult(
        params=params,
        const=const_vec,
        errors=errors,
        sigma=Sigma,
        r_squared=r2,
        vcv=vcv,
        param_vec=param_vec,
        lags=lag_list,
        n_obs=n,
        n_vars=K,
        log_likelihood=float(ll),
    )

grangercause ¶

grangercause(y: FloatArray, lags: int | list[int] = 1, include_const: bool = True, het: bool = True, uncorr: bool = False, method: Literal['lr', 'lm', 'wald'] = 'lr') -> GCResult

Granger causality testing in a VAR.

stat[i,j] tests H0: lags of y_j do not Granger-cause y_i.

Parameters:

Name	Type	Description	Default
`y`	`same as vectorar`		required
`lags`	`same as vectorar`		required
`include_const`	`same as vectorar`		required
`het`	`same as vectorar`		required
`uncorr`	`same as vectorar`		required
`method`	`'lr' \| 'lm' \| 'wald'`	"lr" — likelihood ratio (Chi2, robust if het=True) "lm" — score/LM test "wald" — Wald test using VCV from vectorar	`'lr'`

Returns:

Type	Description
`GCResult with (K,K) matrices of statistics and p-values.`

Source code in src/mfe/timeseries/var.py

def grangercause(
    y: FloatArray,
    lags: int | list[int] = 1,
    include_const: bool = True,
    het: bool = True,
    uncorr: bool = False,
    method: Literal["lr", "lm", "wald"] = "lr",
) -> GCResult:
    """
    Granger causality testing in a VAR.

    stat[i,j] tests H0: lags of y_j do not Granger-cause y_i.

    Parameters
    ----------
    y, lags, include_const, het, uncorr : same as vectorar
    method : "lr" | "lm" | "wald"
        "lr"   — likelihood ratio (Chi2, robust if het=True)
        "lm"   — score/LM test
        "wald" — Wald test using VCV from vectorar

    Returns
    -------
    GCResult with (K,K) matrices of statistics and p-values.
    """
    Y = np.asarray(y, dtype=np.float64)
    T, K = Y.shape

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

    P = len(lag_list)
    df = P  # number of restrictions per (i,j) pair

    res_unr = vectorar(Y, lags=lag_list, include_const=include_const, het=het, uncorr=uncorr)

    stat_mat = np.full((K, K), np.nan, dtype=np.float64)
    pval_mat = np.full((K, K), np.nan, dtype=np.float64)

    Y_adj, X_full, max_lag = _build_regressor_matrix(Y, lag_list, include_const)
    n = Y_adj.shape[0]
    n_regs = X_full.shape[1]
    n_const = 1 if include_const else 0

    for i in range(K):       # caused variable (equation)
        for j in range(K):   # causing variable (excluded)
            if i == j:
                stat_mat[i, j] = np.nan
                pval_mat[i, j] = np.nan
                continue

            # Build restricted X: drop columns for lags of y_j
            # Column layout: [const?] [lag1: K cols] [lag2: K cols] ...
            drop_cols = []
            for p in range(P):
                col_start = n_const + p * K
                drop_cols.append(col_start + j)   # j-th variable in lag p block
            keep_cols = [c for c in range(n_regs) if c not in drop_cols]
            X_r = X_full[:, keep_cols]

            y_i = Y_adj[:, i]

            if method == "wald":
                # Wald: R beta = 0 using VCV from unrestricted
                # Extract VCV block for equation i
                eqs_per_col = n_regs
                vcv_i = res_unr.vcv[i*eqs_per_col:(i+1)*eqs_per_col,
                                     i*eqs_per_col:(i+1)*eqs_per_col]
                # Restriction matrix R: picks rows corresponding to j's lags
                R = np.zeros((df, n_regs), dtype=np.float64)
                for p_idx, col in enumerate(drop_cols):
                    R[p_idx, col] = 1.0
                # Param vector for equation i
                B_i = res_unr.param_vec[i*n_regs:(i+1)*n_regs]
                Rb = R @ B_i
                try:
                    W = float(n * Rb @ np.linalg.solve(R @ vcv_i @ R.T * n, Rb))
                except np.linalg.LinAlgError:
                    W = np.nan
                stat_mat[i, j] = W
                pval_mat[i, j] = float(1 - stats.chi2.cdf(W, df=df)) if np.isfinite(W) else np.nan

            else:
                # LR or LM: compare restricted and unrestricted models for equation i
                # Unrestricted
                XTX_inv_u = np.linalg.inv(X_full.T @ X_full)
                b_u = XTX_inv_u @ (X_full.T @ y_i)
                e_u = y_i - X_full @ b_u
                s2_u = float(e_u @ e_u) / n

                # Restricted
                try:
                    XTX_inv_r = np.linalg.inv(X_r.T @ X_r)
                except np.linalg.LinAlgError:
                    continue
                b_r = XTX_inv_r @ (X_r.T @ y_i)
                e_r = y_i - X_r @ b_r
                s2_r = float(e_r @ e_r) / n

                if method == "lr":
                    if het:
                        # Robust LR-class: based on scores under null, VCV under alt
                        scores = X_full[:, drop_cols] * e_u[:, None]
                        B_meat = scores.T @ scores / n
                        B_bread = X_full[:, drop_cols].T @ X_full[:, drop_cols] / n
                        try:
                            B_inv = np.linalg.inv(B_bread)
                            S_inv = np.linalg.inv(B_inv @ B_meat @ B_inv / n)
                        except np.linalg.LinAlgError:
                            continue
                        s_bar = scores.mean(axis=0)
                        LR = float(n * s_bar @ S_inv @ s_bar)
                    else:
                        LR = float(n * (np.log(s2_r) - np.log(s2_u)))
                    stat_mat[i, j] = LR
                    pval_mat[i, j] = float(1 - stats.chi2.cdf(LR, df=df))

                else:  # lm
                    # LM: regress e_r on X_full, R^2 * n
                    b_aux = np.linalg.lstsq(X_full, e_r, rcond=None)[0]
                    e_aux = e_r - X_full @ b_aux
                    ss_res_aux = float(e_aux @ e_aux)
                    ss_tot_aux = float(e_r @ e_r)
                    r2_aux = 1.0 - ss_res_aux / ss_tot_aux if ss_tot_aux > 0 else 0.0
                    LM = float(n * r2_aux)
                    stat_mat[i, j] = LM
                    pval_mat[i, j] = float(1 - stats.chi2.cdf(LM, df=df))

    return GCResult(
        statistics=stat_mat,
        p_values=pval_mat,
        method=method,
        n_obs=n,
        n_vars=K,
    )

impulse_response ¶

impulse_response(y: FloatArray, lags: int | list[int] = 1, horizon: int = 12, decomp: Literal['unit', 'cholesky', 'spectral'] = 'cholesky', include_const: bool = True, het: bool = True, uncorr: bool = False) -> IRFResult

Impulse response functions for a VAR(P) with delta-method standard errors.

Parameters:

Name	Type	Description	Default
`y`	`FloatArray`		required
`lags`	`int \| list[int]`		`1`
`horizon`	`int`		`12`
`decomp`	`Literal['unit', 'cholesky', 'spectral']`	"unit" — unit shocks (unscaled), i.e. P0 = I_K "cholesky" — Cholesky of Sigma (lower triangular), recursive identification "spectral" — symmetric square root of Sigma (spectral decomposition)	`'cholesky'`
`het`	`VCV options for standard error computation`		`True`
`uncorr`	`VCV options for standard error computation`		`True`

Returns:

Type	Description
`IRFResult`	.responses : (K, K, H+1) — responses[response_var, shock_var, h] .std_errors : (K, K, H+1) delta-method std errors

Source code in src/mfe/timeseries/var.py

def impulse_response(
    y: FloatArray,
    lags: int | list[int] = 1,
    horizon: int = 12,
    decomp: Literal["unit", "cholesky", "spectral"] = "cholesky",
    include_const: bool = True,
    het: bool = True,
    uncorr: bool = False,
) -> IRFResult:
    """
    Impulse response functions for a VAR(P) with delta-method standard errors.

    Parameters
    ----------
    y        : (T, K) data
    lags     : VAR lag order or list
    horizon  : number of periods H; returns H+1 responses (including period 0)
    decomp   : shock decomposition
        "unit"      — unit shocks (unscaled), i.e. P0 = I_K
        "cholesky"  — Cholesky of Sigma (lower triangular), recursive identification
        "spectral"  — symmetric square root of Sigma (spectral decomposition)
    het, uncorr : VCV options for standard error computation

    Returns
    -------
    IRFResult
        .responses  : (K, K, H+1) — responses[response_var, shock_var, h]
        .std_errors : (K, K, H+1) delta-method std errors
    """
    Y = np.asarray(y, dtype=np.float64)
    T, K = Y.shape

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

    res = vectorar(Y, lags=lag_list, include_const=include_const, het=het, uncorr=uncorr)
    P = len(lag_list)
    Sigma = res.sigma

    # Shock decomposition matrix P0
    if decomp == "unit":
        P0 = np.eye(K)
    elif decomp == "cholesky":
        try:
            P0 = np.linalg.cholesky(Sigma)    # lower triangular
        except np.linalg.LinAlgError:
            warnings.warn("Cholesky failed; using spectral decomposition.", RuntimeWarning)
            P0 = _spectral_sqrt(Sigma)
    elif decomp == "spectral":
        P0 = _spectral_sqrt(Sigma)
    else:
        raise ValueError(f"decomp must be 'unit', 'cholesky', or 'spectral', got '{decomp}'")

    # Companion form: convert VAR(P) to VAR(1) state A of size (P*K, P*K)
    # A = [[Phi_1, Phi_2, ..., Phi_P],
    #      [I_K,   0,    ..., 0     ],
    #      [0,     I_K,  ..., 0     ],
    #      ...                       ]
    max_lag = max(lag_list)
    A = np.zeros((max_lag * K, max_lag * K), dtype=np.float64)

    for p_idx, p in enumerate(lag_list):
        # res.params[p_idx] is (K, K) Phi_p
        A[:K, p_idx * K: (p_idx + 1) * K] = res.params[p_idx]
    for block in range(1, max_lag):
        A[block * K: (block + 1) * K, (block - 1) * K: block * K] = np.eye(K)

    # IRF via MA(inf) representation: Psi_h = J A^h J' where J = [I_K | 0]
    J = np.zeros((K, max_lag * K), dtype=np.float64)
    J[:K, :K] = np.eye(K)

    responses  = np.empty((K, K, horizon + 1), dtype=np.float64)
    std_errors = np.empty((K, K, horizon + 1), dtype=np.float64)

    Ah = np.eye(max_lag * K, dtype=np.float64)
    for h in range(horizon + 1):
        Psi_h = J @ Ah @ J.T        # (K, K)
        responses[:, :, h] = Psi_h @ P0
        Ah = Ah @ A

    # Delta-method std errors via asymptotic approximation
    # For each horizon h, Var(vec(Psi_h * P0)) via chain rule on A^h
    # We use numerical differentiation for generality
    std_errors = _irf_std_errors(res, lag_list, horizon, P0, J, A, K, max_lag)

    return IRFResult(
        responses=responses,
        std_errors=std_errors,
        lags=len(lag_list),
        horizon=horizon,
        decomp=decomp,
    )

mfe.timeseries¶

mfe.timeseries ¶

VARResult dataclass ¶

GCResult dataclass ¶

IRFResult dataclass ¶

BNResult dataclass ¶

vectorar ¶

grangercause ¶

impulse_response ¶

beveridge_nelson ¶

BNResult dataclass ¶

beveridge_nelson ¶

var ¶

VARResult dataclass ¶

GCResult dataclass ¶

IRFResult dataclass ¶

vectorar ¶

grangercause ¶

impulse_response ¶

VARResult `dataclass` ¶

GCResult `dataclass` ¶

IRFResult `dataclass` ¶

BNResult `dataclass` ¶

BNResult `dataclass` ¶

VARResult `dataclass` ¶

GCResult `dataclass` ¶

IRFResult `dataclass` ¶