Econometrics Toolbox: Methods guide
← Overview · Getting started · Modeling guide · Methods guide · Examples · Reference
Mathematical foundations
The Econometrics Toolbox implements rigorous statistical algorithms for time series estimation and inference. This guide documents the mathematical formulations, state-space models, parameter estimators, and diagnostic testing routines implemented across the workbench.
1. ARIMA & SARIMAX estimation
Univariate ARIMA\((p,d,q)\) specification
A non-seasonal autoregressive integrated moving average model for a scalar series \(Y_t\) is defined as:
\[\phi(L) (1 - L)^d Y_t = c + \theta(L) \varepsilon_t\]
where: - \(L\) is the lag operator (\(L^k Y_t = Y_{t-k}\)); - \((1 - L)^d = \Delta^d\) is the difference operator of order \(d\); - \(\phi(L) = 1 - \phi_1 L - \phi_2 L^2 - \dots - \phi_p L^p\) is the autoregressive polynomial; - \(\theta(L) = 1 + \theta_1 L + \theta_2 L^2 + \dots + \theta_q L^q\) is the moving average polynomial; - \(c\) is a drift parameter; - \(\varepsilon_t \overset{iid}{\sim} \mathcal{N}(0, \sigma^2)\) represents Gaussian white-noise innovations.
Seasonal multiplicative SARIMAX\((p,d,q)\times(P,D,Q)_s\) with exogenous regressors
When seasonality and exogenous indicators (\(X_t\)) are present, the general multiplicative formulation is:
\[\phi(L)\Phi(L^s) (1 - L)^d (1 - L^s)^D (Y_t - X_t \beta) = \theta(L)\Theta(L^s) \varepsilon_t\]
where: - \(\Phi(L^s) = 1 - \Phi_1 L^s - \dots - \Phi_P L^{Ps}\) is the seasonal autoregressive polynomial; - \(\Theta(L^s) = 1 + \Theta_1 L^s + \dots + \Theta_Q L^{Qs}\) is the seasonal moving average polynomial; - \(s\) is the seasonal period (\(s=4\) for quarterly, \(s=12\) for monthly); - \(D\) is the seasonal differencing order; - \(X_t\) is a row vector of exogenous regressors and \(\beta\) is the corresponding coefficient vector.
State-space representation & Kalman filter
The notebook casts SARIMAX models into state-space form:
\[\begin{aligned} \alpha_{t+1} &= T \alpha_t + R \eta_t, \quad &\eta_t \sim \mathcal{N}(0, Q) \quad &\text{(State transition equation)} \\ Y_t &= Z \alpha_t + d_t + \varepsilon_t, \quad &\varepsilon_t \sim \mathcal{N}(0, H) \quad &\text{(Observation equation)} \end{aligned}\]
The Kalman filter recursively updates state predictions (\(\hat{\alpha}_{t|t-1}\)), error covariances (\(P_{t|t-1}\)), and forecast errors (\(v_t = Y_t - Z \hat{\alpha}_{t|t-1} - d_t\)) with variance \(F_t = Z P_{t|t-1} Z' + H\).
Exact Gaussian maximum likelihood
Rather than approximating the initial conditions via conditional sum-of-squares (CSS), the exact Gaussian log-likelihood is maximized over parameters \(\theta = (\phi, \theta, \Phi, \Theta, \beta, \sigma^2)\):
\[\ln L(\theta) = -\frac{N}{2} \ln(2\pi) - \frac{1}{2} \sum_{t=1}^N \left( \ln |F_t| + v_t' F_t^{-1} v_t \right)\]
Parameter optimization is executed using quasi-Newton algorithms (L-BFGS). The parameter covariance matrix is evaluated from the inverted Hessian:
\[\widehat{\text{Var}}(\hat{\theta}) = \left[ -\frac{\partial^2 \ln L}{\partial \theta \partial \theta'} \right]^{-1}\]
yielding asymptotic standard errors \(SE(\hat{\theta}_j) = \sqrt{\widehat{\text{Var}}(\hat{\theta})_{jj}}\) and \(z\)-statistics \(z = \hat{\theta}_j / SE(\hat{\theta}_j)\).
2. Vector Autoregression (VAR)
Specification
A \(K\)-variable Vector Autoregression of order \(p\), denoted \(\text{VAR}(p)\), models the joint dynamic evolution of an endogenous vector \(Y_t = [Y_{1,t}, Y_{2,t}, \dots, Y_{K,t}]'\):
\[Y_t = c + A_1 Y_{t-1} + A_2 Y_{t-2} + \dots + A_p Y_{t-p} + u_t\]
where: - \(c\) is a \(K \times 1\) vector of intercepts; - \(A_i\) are \(K \times K\) coefficient matrices; - \(u_t\) is a \(K \times 1\) vector of unobservable zero-mean white noise with covariance matrix \(\Sigma_u = \mathbb{E}[u_t u_t']\).
In the supplied treasury_rates_var fixture, \(K=2\) (\(Y_{1,t} = \text{ShortRate}\), \(Y_{2,t} = \text{LongRate}\)) and \(p=1\):
\[\begin{bmatrix} Y_{1,t} \\ Y_{2,t} \end{bmatrix} = \begin{bmatrix} c_1 \\ c_2 \end{bmatrix} + \begin{bmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{bmatrix} \begin{bmatrix} Y_{1,t-1} \\ Y_{2,t-1} \end{bmatrix} + \begin{bmatrix} u_{1,t} \\ u_{2,t} \end{bmatrix}\]
Estimation & stability
Each equation is estimated via Ordinary Least Squares (OLS), which is numerically equivalent to Generalized Least Squares (GLS) when regressors are identical across equations.
The system is covariance-stationary if all eigenvalues \(\lambda\) of the companion matrix \(F\) lie strictly inside the complex unit circle (\(|\lambda| < 1\)):
\[F = \begin{bmatrix} A_1 & A_2 & \dots & A_{p-1} & A_p \\ I_K & 0 & \dots & 0 & 0 \\ 0 & I_K & \dots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & I_K & 0 \end{bmatrix}\]
Granger causality test
To test whether variable \(Y_2\) Granger-causes variable \(Y_1\), the null hypothesis \(H_0: A_{12,1} = A_{12,2} = \dots = A_{12,p} = 0\) is tested using a Wald \(F\)-statistic:
\[F = \frac{(SSR_R - SSR_U) / p}{SSR_U / (N - Kp - 1)}\]
where \(SSR_R\) is the restricted sum of squared residuals (excluding \(Y_2\) lags) and \(SSR_U\) is the unrestricted sum of squared residuals.
3. Cointegration & Error Correction (ECM)
The cointegration concept
Two non-stationary series \(Y_t \sim I(1)\) and \(X_t \sim I(1)\) are cointegrated of order \((1,1)\), denoted \(CI(1,1)\), if there exists a linear combination:
\[e_t = Y_t - \beta_0 - \beta_1 X_t\]
such that \(e_t \sim I(0)\) is a stationary process. The vector \([1, -\beta_1]'\) represents the long-run equilibrium relationship.
Engle-Granger two-step procedure
- Step 1 (Static OLS Regression): Estimate the cointegrating equation via OLS: \[Y_t = \beta_0 + \beta_1 X_t + e_t\] and obtain residuals \(\hat{e}_t = Y_t - \hat{\beta}_0 - \hat{\beta}_1 X_t\).
- Step 2 (Residual Unit Root Test): Run an Augmented Dickey-Fuller test on \(\hat{e}_t\): \[\Delta \hat{e}_t = \gamma \hat{e}_{t-1} + \sum_{j=1}^k \psi_j \Delta \hat{e}_{t-j} + v_t\] Because \(\hat{e}_t\) is an estimated rather than observed series, standard Dickey-Fuller critical values do not apply; the notebook computes MacKinnon cointegration \(p\)-values via
statsmodels.tsa.stattools.coint.
Granger Representation Theorem & Dynamic ECM
By the Granger Representation Theorem, if \(Y_t\) and \(X_t\) are cointegrated, they can be represented as an Error Correction Model (ECM):
\[\Delta Y_t = c_0 + \alpha \hat{e}_{t-1} + \sum_{i=1}^p \delta_i \Delta Y_{t-i} + \sum_{j=0}^q \theta_j \Delta X_{t-j} + \varepsilon_t\]
where: - \(\hat{e}_{t-1} = Y_{t-1} - \hat{\beta}_0 - \hat{\beta}_1 X_{t-1}\) is the error correction term; - \(\alpha\) is the speed of adjustment coefficient. Stability requires \(\alpha < 0\): when \(Y_{t-1}\) exceeds its long-run equilibrium (\(\hat{e}_{t-1} > 0\)), \(\Delta Y_t\) contracts in the next period to restore equilibrium.
In the supplied consumption_income_coint fixture, \(\alpha = -0.8051\) (\(p < 0.05\)), indicating rapid return to long-run equilibrium following income shocks.
4. Stationarity & unit root testing
Difference-stationary vs. trend-stationary
- Difference-Stationary (\(I(1)\)): Contains a stochastic trend with infinite memory. Shocks permanently alter the level: \[Y_t = \mu + Y_{t-1} + \varepsilon_t \implies Y_t = Y_0 + \mu t + \sum_{i=1}^t \varepsilon_i\]
- Trend-Stationary (\(I(0)\) with drift): Fluctuates around a deterministic trend line with finite memory: \[Y_t = \alpha + \beta t + \varepsilon_t\]
Augmented Dickey-Fuller (ADF) test
The ADF test evaluates the null hypothesis of a unit root (\(H_0: \gamma = 0\)) against the stationary alternative (\(H_1: \gamma < 0\)):
\[\Delta Y_t = \alpha + \beta t + \gamma Y_{t-1} + \sum_{j=1}^k \delta_j \Delta Y_{t-j} + u_t\]
The test statistic is the \(t\)-ratio:
\[\tau = \frac{\hat{\gamma}}{SE(\hat{\gamma})}\]
Because \(\tau\) does not follow a Student’s \(t\)-distribution under non-stationarity, the notebook evaluates finite-sample MacKinnon approximate \(p\)-values. A \(p\)-value \(< 0.05\) rejects the unit root, confirming stationarity.
KPSS stationarity test
The Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test flips the null and alternative hypotheses: - \(H_0: Y_t\) is level/trend stationary. - \(H_1: Y_t\) has a unit root.
The series is decomposed into \(Y_t = \xi t + r_t + \varepsilon_t\), where \(r_t = r_{t-1} + u_t\) with \(u_t \sim \mathcal{N}(0, \sigma_u^2)\). Under \(H_0\), \(\sigma_u^2 = 0\). The LM test statistic is:
\[\text{KPSS} = \frac{\sum_{t=1}^N S_t^2}{N^2 s^2(l)}\]
where \(S_t = \sum_{i=1}^t e_i\) is the partial sum of residuals from regressing \(Y_t\) on a constant, and \(s^2(l)\) is a Newey-West long-run variance estimator. A \(p\)-value \(\ge 0.05\) fails to reject stationarity.
5. Residual diagnostic battery
Post-estimation residual diagnostics ensure that the fitted model adequately captures all systematic variation in the historical series:
\[\hat{\varepsilon}_t = Y_t - \hat{Y}_t\]
Ljung-Box \(Q\) portmanteau test
Tests whether the first \(m\) autocorrelation coefficients of the residuals are jointly zero (\(H_0: \rho_1 = \dots = \rho_m = 0\)):
\[Q = N(N + 2) \sum_{k=1}^m \frac{\hat{\rho}_k^2}{N - k} \sim \chi^2(m)\]
where \(\hat{\rho}_k\) is the sample autocorrelation of \(\hat{\varepsilon}_t\) at lag \(k\). A \(p\)-value \(\ge 0.05\) confirms the absence of residual autocorrelation.
Jarque-Bera normality test
Tests whether residuals have the skewness (\(S\)) and kurtosis (\(K\)) matching a Gaussian distribution (\(H_0: S = 0 \text{ and } K = 3\)):
\[JB = \frac{N}{6} \left( S^2 + \frac{(K - 3)^2}{4} \right) \sim \chi^2(2)\]
where \(S = \frac{\frac{1}{N}\sum \hat{\varepsilon}_t^3}{(\frac{1}{N}\sum \hat{\varepsilon}_t^2)^{3/2}}\) and \(K = \frac{\frac{1}{N}\sum \hat{\varepsilon}_t^4}{(\frac{1}{N}\sum \hat{\varepsilon}_t^2)^2}\).
Engle’s ARCH LM test
Tests for autoregressive conditional heteroskedasticity (volatility clustering) in residuals by regressing squared residuals on their own lags:
\[\hat{\varepsilon}_t^2 = \alpha_0 + \alpha_1 \hat{\varepsilon}_{t-1}^2 + \dots + \alpha_q \hat{\varepsilon}_{t-q}^2 + v_t\]
The test statistic is:
\[LM = N \cdot R^2 \sim \chi^2(q)\]
Under \(H_0: \alpha_1 = \dots = \alpha_q = 0\), residuals are homoskedastic. A \(p\)-value \(\ge 0.05\) confirms constant error variance.
6. Information criteria
To penalize overfitting when choosing lag orders \(p\) and \(q\), the toolbox reports four criterion metrics:
| Criterion | Formula | Asymptotic Property |
|---|---|---|
| Akaike Information Criterion (AIC) | \(\text{AIC} = -2 \ln L + 2k\) | Efficient in finite samples; tends to slightly overparameterize as \(N \to \infty\) |
| Bayesian Information Criterion (BIC) | \(\text{BIC} = -2 \ln L + k \ln(N)\) | Strongly consistent: selects the true lag order with probability 1 as \(N \to \infty\) |
| Hannan-Quinn Criterion (HQIC) | \(\text{HQIC} = -2 \ln L + 2k \ln(\ln(N))\) | Consistent; intermediate penalty between AIC and BIC |
| Log-Likelihood (\(\ln L\)) | \(\ln L(\hat{\theta})\) | Unpenalized goodness-of-fit |
where \(k\) is the number of estimated parameters and \(N\) is the sample size.
MATLAB Econometrics Toolbox capability comparison
The table below contrasts the capabilities of the Boardflare Econometrics Toolbox against the commercial MATLAB Econometrics Toolbox:
| Capability Area | MATLAB Function | Boardflare Econometrics Toolbox | Implementation Notes |
|---|---|---|---|
| ARIMA Estimation | arima, estimate |
statsmodels.tsa.arima.model.ARIMA |
Both employ exact state-space Kalman filter MLE. Boardflare runs directly inside Excel without client licenses. |
| Seasonal SARIMAX | arima with SAR / SMA / X |
statsmodels.tsa.statespace.sarimax.SARIMAX |
Multiplicative seasonal orders and exogenous indicators fully supported in both. |
| Vector Autoregression | varm, estimate |
statsmodels.tsa.api.VAR |
Multi-equation OLS estimation, companion matrix stability, and multi-step joint forecasting. |
| Granger Causality | Custom Wald test script | VARResults.test_causality |
Automated Wald \(F\)-test reporting test statistics and \(p\)-values. |
| Cointegration Testing | egcitest |
statsmodels.tsa.stattools.coint |
Engle-Granger two-step procedure with MacKinnon asymptotic \(p\)-values. |
| Error Correction Models | vecm |
Custom two-step OLS dynamic ECM | Recovers long-run cointegrating vector and dynamic adjustment speed \(\alpha\). |
| Unit Root ADF | adftest |
statsmodels.tsa.stattools.adfuller |
Dickey-Fuller regression with configurable maximum lag and AIC lag selection. |
| Stationarity KPSS | kpsstest |
statsmodels.tsa.stattools.kpss |
Tests null of stationarity using Newey-West long-run variance estimation. |
| Autocorrelation | lbqtest |
statsmodels.stats.diagnostic.acorr_ljungbox |
Portmanteau \(Q\)-statistic against \(\chi^2\) distribution. |
| Normality | jbtest |
statsmodels.stats.stattools.jarque_bera |
Skewness-kurtosis joint test with \(\chi^2(2)\) distribution. |
| Heteroskedasticity | archtest |
statsmodels.stats.diagnostic.het_arch |
Engle’s LM test for conditional autoregressive variance. |
| Information Criteria | aicbic |
Built-in AIC, BIC, HQIC | Automated computation across all estimated model families. |
| GARCH Volatility | garch, egarch |
Planned (Phase E / Future) | MATLAB provides built-in GARCH; Boardflare plan targets arch library integration. |
| State-Space Filtering | ssm |
Exploring | State-space unobserved components models exploring integration. |