ngme VAR(1) bivariate model specification (Cayley re-parameterization)

Builds a Vector Autoregressive order-1 (VAR(1)) operator for bivariate time-series data. The model follows: $$Y_t = A\,Y_{t-1} + \varepsilon_t, \quad A = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{pmatrix}$$

Usage

var1(mesh = NULL, p1 = 0, p2 = 1, p3 = 0, p4 = 1)

Arguments

mesh: Integer vector or inla.mesh.1d object for the time mesh (length $T$). Pass NULL for deferred construction.
p1: Unconstrained skew-symmetric parameter ($\in \mathbb{R}$); controls the rotation part of the Cayley mapping. Default 0.
p2, p3, p4: Unconstrained lower-triangular Cholesky parameters ($\in \mathbb{R}$); control the contraction (positive-definite) part. Defaults p2 = 1, p3 = 0, p4 = 1.

Value

An ngme_operator (or ngme_operator_def when mesh = NULL) suitable for use inside f().

Cayley re-parameterization

To guarantee stationarity ($\rho(A) < 1$) at every SGD step, the $2 \times 2$ coefficient matrix $A$ is obtained via the Cayley transform: $$A = (I + S)(I - S)^{-1}$$ where $S$ is constructed from four unconstrained real parameters $(p_1, p_2, p_3, p_4) \in \mathbb{R}^4$ as $S = J - R$:

Skew-symmetric part (rotation):: $J = \begin{pmatrix} 0 & p_1 \\ -p_1 & 0 \end{pmatrix}$
Positive-definite part (contraction):: $R = L L^{\top} + \varepsilon I$, where $L = \begin{pmatrix} p_2 & 0 \\ p_3 & p_4 \end{pmatrix}$ and $\varepsilon = 10^{-5}$.

Because $S + S^{\top} = -2R$ is negative definite, all eigenvalues of $S$ have strictly negative real parts, and the Cayley transform then guarantees $|\lambda_i(A)| < 1$ for every $i$.

The default values $(p_1, p_2, p_3, p_4) = (0, 1, 0, 1)$ give $A \approx 0$ (no auto-regression at initialization).

Precision operator

The 2T x 2T precision operator is: $$K = M_0 + a_{11}\,M_{11} + a_{22}\,M_{22} + a_{12}\,M_{12} + a_{21}\,M_{21}$$ where the $M_{ij}$ are fixed sparse block matrices constructed from the $T \times T$ first-order difference matrix $C_T$.

Examples

if (FALSE) { # \dontrun{
n_obs <- 200
dat <- data.frame(
  y     = c(y1_series, y2_series),
  idx   = rep(seq_len(n_obs), 2),
  group = factor(rep(c("y1", "y2"), each = n_obs))
)

fit <- ngme(
  y ~ 0 + f(idx, model = var1(mesh = 1:n_obs), group = group,
            noise = noise_nig()),
  data    = dat,
  family  = noise_normal(sigma = 0.01, fix_sigma = TRUE),
  control_opt = control_opt(iterations = 1000, n_parallel_chain = 4)
)
print(fit)
} # }