---
title: "Systems of predictive regressions"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Systems of predictive regressions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

```{r setup, message = FALSE}
library(ivx)
```

`ivx_sys()` estimates several predictive regressions jointly — for instance
the returns on a set of pricing factors on the same predictors, the empirical
setting of Kostakis, Magdalinos & Stamatogiannis (2023) — and tests
restrictions across equations with a single IVX-Wald statistic.

## Model and statistic

$$
y_t = \mu + A\, x_{t-1} + \varepsilon_t, \qquad y_t \in \mathbb{R}^m, \; x_t \in \mathbb{R}^r,
$$

with the predictors as in `vignette("ivx")` and $\varepsilon_t$ a vector
martingale difference with covariance $\Sigma$. The IVX estimator of the
$m \times r$ matrix $A$ (KMS 2023, eq. 15) is

$$
\tilde A_K = Y(K)' Z \,\big[ X(K)' Z \big]^{-1},
$$

and the covariance of $\mathrm{vec}(\tilde A_K)$ has the Kronecker form of
eq. (23),

$$
\tilde Q_K = \big[ (Z'X(K))^{-1} \otimes I_m \big]\, M_K\,
             \big[ (X(K)'Z)^{-1} \otimes I_m \big], \qquad
M_K = Z(K)'Z(K) \otimes \hat\Sigma \;-\; n_K\, \bar z(K)\bar z(K)' \otimes \hat\Sigma_{FM},
$$

where $\hat\Sigma$ is the covariance of the one-step OLS residuals and
$\hat\Sigma_{FM} = \hat\Sigma - \hat\Omega_{\varepsilon v}' \hat\Omega_{vv}^{-1} \hat\Omega_{\varepsilon v}$
its long-run-corrected analogue. For $H \mathrm{vec}(A) = h$,

$$
W = \big(H\mathrm{vec}(\tilde A_K) - h\big)' \big[H \tilde Q_K H'\big]^{-1} \big(H\mathrm{vec}(\tilde A_K) - h\big)
\;\to\; \chi^2_{\mathrm{rank}(H)} .
$$

`ivx_sys()` reports the joint test of $A = 0$ ($mr$ degrees of freedom), one
test per equation ($r$ degrees of freedom) and one per coefficient. For
$m = 1$ every number equals the output of `ivx()`, at any horizon.

```{r}
s <- ivx_sys(cbind(Ret, DE) ~ DP + TBL, data = kms)
s
summary(s)
```

Long horizons work as for `ivx()`:

```{r}
s12 <- ivx_sys(cbind(Ret, DE) ~ DP + TBL, data = kms, horizon = 12)
s12$Wald_Eq
```

`vcov()` returns the $mr \times mr$ matrix of $\mathrm{vec}(A)$, column-major
(all responses for the first predictor, then the second, ...), with names
`response:predictor`; custom restrictions can be tested from it directly.

```{r}
V <- vcov(s)
a <- as.vector(coef(s))
names(a) <- rownames(V)
# is the effect of DP the same in both equations?
H <- matrix(c(1, -1, 0, 0), 1)
drop(t(H %*% a) %*% solve(H %*% V %*% t(H)) %*% (H %*% a))
```

## Validation

* For one response, `ivx_sys()` reproduces `ivx()` exactly (tested for
  horizons 1 and 4), and the first equation's Wald statistic in a system
  equals the univariate joint Wald statistic.
* Monte Carlo with two responses, a unit-root predictor and endogeneity
  correlations $-0.95$ and $0.5$, $n = 500$: empirical size of the joint 5%
  test 4.4%, of the equation tests 5.6% and 5.0%.

## Caveats

* The example uses `DE` as a second response only to illustrate the syntax;
  the theory assumes stationary responses.
* All equations share the same predictors and the same instrument; equation
  specific regressors are not supported.
* No Eicker–White option, no bootstrap and no `ivx_ar`/`ivx_ra` counterpart
  for systems.

## References

* Kostakis, A., Magdalinos, T., & Stamatogiannis, M. P. (2023). Taking stock
  of long-horizon predictability tests: Are factor returns predictable?
  *Journal of Econometrics*, 237(2), 105380.
* Magdalinos, T. (2022). Least squares and IVX limit theory in systems of
  predictive regressions with GARCH innovations. *Econometric Theory*, 38(5),
  875–912.