From fe71033e6f442596835c9a99b0a77b3ee2f2defd Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 09:26:10 +0000 Subject: [PATCH 01/16] Add implementation plan and prototype for complex generalised normal (dcgnorm) Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- .Rbuildignore | 1 + dev/PLAN-dcgnorm.md | 122 ++++++++++++++++++++++++++++++++++++++++ dev/prototype-dcgnorm.R | 37 ++++++++++++ 3 files changed, 160 insertions(+) create mode 100644 dev/PLAN-dcgnorm.md create mode 100644 dev/prototype-dcgnorm.R diff --git a/.Rbuildignore b/.Rbuildignore index 9fb1e4e..9048bdd 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -9,3 +9,4 @@ man-roxygen ^revdep$ cran-comments.md ^CRAN-SUBMISSION$ +^dev$ diff --git a/dev/PLAN-dcgnorm.md b/dev/PLAN-dcgnorm.md new file mode 100644 index 0000000..49b5ade --- /dev/null +++ b/dev/PLAN-dcgnorm.md @@ -0,0 +1,122 @@ +# Plan: complex generalised normal distribution (`dcgnorm`) in `complex` + +Branch: `dcgnorm` (from `master` at `a8b660c`, v1.0.2). Prototype with numerical checks: +`dev/prototype-dcgnorm.R`. `dev/` is in `.Rbuildignore`. + +Motivation: residuals of wind data have much heavier tails than normal (see the cARIMA wind +project). The complex normal likelihood in `clm()` is the only likelihood available. Goal: add the +complex generalised normal (CGN) with shape β estimated, and only ONE extra estimated parameter, +by concentrating the scale and pseudo-scale out analytically. + +## 1. The distribution + +Parameters: location μ (complex), scale α² > 0 (real), pseudo-scale ψ² (complex, |ψ²| < α²), +shape β > 0. For e = z − μ: + + Q(e) = 2 (α²|e|² − Re(Conj(ψ²) e²)) / (α⁴ − |ψ²|²) + f(z) = 2β / (π Γ(1/β) 2^(1/β) sqrt(α⁴ − |ψ²|²)) · exp(−½ Q^β) + +- β = 1: complex normal with σ² = α², ς² = ψ² (checked in the prototype: max log-density difference from `dcnorm` is 2e-16). +- β < 1 heavier tails (β = 0.5 Laplace-like), β → ∞ uniform on an ellipse. +- Moments: variance σ² = k(β) α², pseudo-variance ς² = k(β) ψ², where + k(β) = 2^(1/β) Γ(2/β) / (2 Γ(1/β)). Hence ς²/σ² = ψ²/α² for all β. +- Radial law: Q^β ~ Gamma(shape = 1/β, scale = 2). Used for simulation and diagnostics. +- Equivalent to the bivariate exponential power distribution (Gómez et al., 1998) and to the + complex generalized Gaussian of Novey, Adalı & Roy (2010, IEEE TSP 58(3)) under the + `cnorm.R` mapping Σ = ½[[α²+Re ψ², Im ψ²],[Im ψ², α²−Re ψ²]]. + +**Parameterisation decision (to confirm):** user-facing functions take `sigma2`, `varsigma2` +(variance and pseudo-variance, as in `dcnorm`) plus `shape`, and convert internally with k(β). +This keeps `sigma.clm()` meaningful and `dcgnorm(shape=1)` == `dcnorm`. Alternative: expose +`scale2`/`pseudoscale2` directly. Recommendation: variance parameterisation. + +## 2. Concentrated likelihood (closed form) + +Write ψ² = ρ α², ρ the complex circularity coefficient (|ρ| < 1). + +1. ρ̂ = mean(e²) / mean(|e|²). Moment estimator, consistent for every β, exact MLE when β = 1. +2. Given ρ̂, qₜ = 2(|eₜ|² − Re(Conj(ρ̂) eₜ²)) / (1 − |ρ̂|²), and the exact MLE of the scale is + α̂² = [ (β / (2T)) Σ qₜ^β ]^(1/β), ψ̂² = ρ̂ α̂². +3. Plugging in gives, since Σ (qₜ/α̂²)^β = 2T/β, + ℓ(B, β) = T [ log(2β) − log π − lgamma(1/β) − log(2)/β − log α̂² − ½ log(1 − |ρ̂|²) − 1/β ]. + With β = 1 this equals the current `clm` cost, T(log 2π + 1 + ½ log det Σ̂). + +Prototype check (T = 20,000, β = 0.5, α² = 2, ψ² = 0.6+0.8i): closed form ℓ = −101981.6 vs +full numeric MLE −101981.5; the β that maximises the closed-form likelihood is 0.50. + +Optional exact-MLE refinement (argument `refine`, default 0): n fixed-point steps +wₜ = Qₜ^(β−1), α² = (β/T) Σ wₜ|eₜ|², ψ² = (β/T) Σ wₜ eₜ², then recompute ℓ with the full +density. Still closed-form arithmetic; the optimiser sees only β. Convergence is proven for +β ≤ 1 (Pascal et al., 2013); for β > 1 cap the iterations and fall back to step 2. + +## 3. Implementation phases + +Each phase ends with `R CMD build` + `R CMD check --no-manual` clean and its tests passing. +Commit per phase. + +### Phase 1 — distribution functions (`R/cgnorm.R`, new) +- `dcgnorm(q, mu=0, sigma2=1, varsigma2=0, shape=1, log=FALSE)`. Default shape 1 = normal. + In greybox `dgnorm` shape 2 = normal (exponent on |x|, not on the quadratic form); document + that shape here equals greybox shape / 2. +- `rcgnorm(n, mu, sigma2, varsigma2, shape)` via the Gamma radial law and a uniform angle. +- Internal helper `cgnormScale(sigma2, varsigma2, shape)` returning (α², ψ²), and its inverse. +- Validate: sigma2 > 0, |varsigma2| < sigma2, shape > 0. +- No `p`/`q` functions in this release (no closed form; `qcnorm` semantics are unusual anyway). +- Roxygen docs in the `cnormal` style; export; NAMESPACE via roxygen. + +### Phase 2 — `clm()` estimation (`R/clm.R`) +- New argument `distribution = c("dcnorm", "dcgnorm")`, used only when `loss = "likelihood"`. + Default keeps current behaviour exactly. +- Shape: pass via `...` as `shape` (NULL = estimate, number = fix), mirroring greybox `alm()` + handling of `other`. Returned in `object$other$shape`. +- `fitter()` (≈ lines 195–247): for `dcgnorm`, compute ρ̂, α̂², ψ̂² as in section 2 and return + `scale` as the 2×2 covariance (via k(β)), so downstream code keeps its meaning. +- Parameter vector: append ONE real parameter log β at the end of `BReal` in `estimator()` + (≈ lines 307–341) and strip it BEFORE the real/imaginary split in `fitter()` (lines 197–201), + which currently assumes every entry is half of a complex coefficient. Update `maxeval`. +- `CF()` (≈ lines 259–305): `dcgnorm` branch = −ℓ(B, β) from section 2. +- Starting values: B from the current `dcnorm` fit logic, log β = 0. +- `nParam` (line 925): add 0.5 when β is estimated (package counts in complex units). +- `logLik` stored as −CF; `AIC`/`AICc` then work via `nparam`/`logLik` methods. + +### Phase 3 — inference (`vcov.clm`, ≈ lines 1041–1215) +- Hessian (line 901) currently over B only. For `dcgnorm` with estimated β compute over (B, log β) + and take the B block of the inverse. The refit call in `vcov.clm` must pass `distribution` + and `shape`. +- Delta method for the SE of β if reported in `summary`. + +### Phase 4 — prediction (`predict.clm`, ≈ lines 1368–1600) +- Point forecasts unchanged. +- Intervals for `dcgnorm`: simulate. Draw ε from `rcgnorm` with fitted (σ², ς², β), push through + the ARIMA recursion for h steps, take quantiles of Re and Im separately (matching the current + per-part intervals). `nsim` argument, default 10000. Return the simulated paths optionally + (needed for energy scores in the wind paper). +- Leave the `dcnorm` interval code as is (only warn as now). + +### Phase 5 — output and diagnostics +- `print.summary.clm`: print distribution and shape (with SE when estimated). +- `plot.clm`: the greybox normal Q–Q plot is wrong for `dcgnorm`; add a Q–Q plot of Qₜ^β against + Gamma(1/β, 2) (a new `which` value), or at least skip the normal Q–Q plot. +- NEWS entry, DESCRIPTION version bump (1.1.0), README mention. + +## 4. Tests (`tests/testthat/`, none exist yet; `testthat` already in Suggests) + +1. `dcgnorm(shape=1)` == `dcnorm` (tolerance 1e-12) for several (sigma2, varsigma2). +2. Density integrates to 1 on a grid for shape ∈ {0.3, 0.5, 1, 2, 5}. +3. `rcgnorm`: sample σ², ς² within 3 SE of targets; KS test of Q^β against Gamma(1/β, 2) p > 0.001 + (fixed seed). +4. Closed-form ℓ equals −Σ log `dcgnorm` evaluated at (α̂², ψ̂²) (identity check of section 2). +5. `clm(y~x, distribution="dcgnorm", shape=1)` reproduces `clm(y~x)` coefficients (1e-6) and logLik. +6. cAR(1) simulated with β = 0.5, T = 2000: β̂ in [0.4, 0.6], AICc(dcgnorm) < AICc(dcnorm). +7. Regression: every existing example in the docs gives unchanged results with the default + `distribution`. + +## 5. Open questions for the author + +1. Variance parameterisation (recommended) or scale parameterisation for `dcgnorm`? +2. `distribution` as a new argument (recommended) or extend `loss` values (e.g. `"likelihood-cgn"`)? +3. Parameter counting: `nparam` counts complex units (one complex coefficient = 1, Σ = 1.5). + AIC/AICc of `clm` therefore penalise half the real-parameter count. Keep the convention (and + document it) or switch to real counts? This affects comparing cARIMA with real VAR by IC. +4. Ship `refine` (exact MLE of ρ) in the first version or later? +5. Name: `dcgnorm` (parallel to greybox `dgnorm`) — confirm. diff --git a/dev/prototype-dcgnorm.R b/dev/prototype-dcgnorm.R new file mode 100644 index 0000000..6309bbb --- /dev/null +++ b/dev/prototype-dcgnorm.R @@ -0,0 +1,37 @@ +suppressMessages(library(complex)) +# complex GG density, scatter (alpha2 real, psi2 complex) +dcgg <- function(z, alpha2, psi2, beta, log=TRUE){ + D <- alpha2^2 - Mod(psi2)^2 + Q <- 2*(alpha2*Mod(z)^2 - Re(Conj(psi2)*z^2))/D + l <- log(2*beta) - log(pi) - lgamma(1/beta) - log(2)/beta - 0.5*log(D) - 0.5*Q^beta + if(log) l else exp(l) +} +rcgg <- function(n, alpha2, psi2, beta){ + S <- matrix(c(alpha2+Re(psi2), Im(psi2), Im(psi2), alpha2-Re(psi2)),2)/2 + L <- t(chol(S)); r <- sqrt(rgamma(n, 1/beta, scale=2)^(1/beta)); ph <- runif(n,0,2*pi) + x <- L %*% rbind(r*cos(ph), r*sin(ph)); complex(real=x[1,], imaginary=x[2,]) +} +# closed-form concentrated logLik: rho from moments, alpha2 exact given rho +concLL <- function(e, beta){ + T <- length(e); rho <- mean(e^2)/mean(Mod(e)^2) + q <- 2*(Mod(e)^2 - Re(Conj(rho)*e^2))/(1-Mod(rho)^2) + a2 <- ((beta/(2*T))*sum(q^beta))^(1/beta) + ll <- T*(log(2*beta) - log(pi) - lgamma(1/beta) - log(2)/beta - log(a2) - 0.5*log(1-Mod(rho)^2) - 1/beta) + list(ll=ll, alpha2=a2, psi2=rho*a2) +} +set.seed(1) +# 1. beta=1 equals dcnorm +z <- rcnorm(5, 0, sigma2=2, varsigma2=0.6+0.8i) +cat("beta=1 vs dcnorm max diff:", max(abs(dcgg(z,2,0.6+0.8i,1) - dcnorm(z,0,2,0.6+0.8i,log=TRUE))), "\n") +# 2. density integrates to 1 (beta=0.5) +g <- seq(-60,60,length.out=1201); zz <- outer(g, 1i*g, "+") +cat("integral beta=0.5:", sum(dcgg(zz,2,0.6+0.8i,0.5,log=FALSE))*diff(g)[1]^2, "\n") +# 3. simulate beta=0.5, compare closed-form vs full numeric MLE of (alpha2, psi2) +e <- rcgg(20000, 2, 0.6+0.8i, 0.5) +cf <- concLL(e, 0.5) +negLL <- function(p) -sum(dcgg(e, exp(p[1]), complex(real=p[2],imaginary=p[3]), 0.5)) +op <- optim(c(log(2),0.5,0.5), negLL, control=list(reltol=1e-12, maxit=5000)) +cat("closed-form: LL", round(cf$ll,2), " alpha2", round(cf$alpha2,4), " psi2", format(round(cf$psi2,4)), "\n") +cat("full MLE : LL", round(-op$value,2), " alpha2", round(exp(op$par[1]),4), " psi2", format(round(complex(real=op$par[2],imaginary=op$par[3]),4)), "\n") +# 4. beta profile from closed form +b <- seq(0.3,1.5,0.01); cat("beta maximising closed-form LL:", b[which.max(sapply(b, function(x) concLL(e,x)$ll))], "(true 0.5)\n") From 1a344dfdb2d7a3ebf5cead0dfe2c0afd37988367 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 09:32:15 +0000 Subject: [PATCH 02/16] Plan: record decisions, add IC parameter-count analysis Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- dev/PLAN-dcgnorm.md | 79 ++++++++++++++++++++++++++++++--------------- 1 file changed, 53 insertions(+), 26 deletions(-) diff --git a/dev/PLAN-dcgnorm.md b/dev/PLAN-dcgnorm.md index 49b5ade..a0ea13b 100644 --- a/dev/PLAN-dcgnorm.md +++ b/dev/PLAN-dcgnorm.md @@ -25,10 +25,11 @@ shape β > 0. For e = z − μ: complex generalized Gaussian of Novey, Adalı & Roy (2010, IEEE TSP 58(3)) under the `cnorm.R` mapping Σ = ½[[α²+Re ψ², Im ψ²],[Im ψ², α²−Re ψ²]]. -**Parameterisation decision (to confirm):** user-facing functions take `sigma2`, `varsigma2` -(variance and pseudo-variance, as in `dcnorm`) plus `shape`, and convert internally with k(β). -This keeps `sigma.clm()` meaningful and `dcgnorm(shape=1)` == `dcnorm`. Alternative: expose -`scale2`/`pseudoscale2` directly. Recommendation: variance parameterisation. +**Parameterisation (DECIDED):** scale parameterisation. Functions take `scale` (α², real) and +`pseudoscale` (ψ², complex) plus `shape` (β). The closed form is simplest in the scale (section 2); +the variance form only adds the factor k(β). At shape = 1, scale = σ² and pseudoscale = ς², so +`dcgnorm(q, mu, scale=s, pseudoscale=p, shape=1)` == `dcnorm(q, mu, sigma2=s, varsigma2=p)`. +`sigma.clm()` keeps returning the covariance of residuals; `object$scale` holds (α², ψ²) for dcgnorm. ## 2. Concentrated likelihood (closed form) @@ -44,10 +45,9 @@ Write ψ² = ρ α², ρ the complex circularity coefficient (|ρ| < 1). Prototype check (T = 20,000, β = 0.5, α² = 2, ψ² = 0.6+0.8i): closed form ℓ = −101981.6 vs full numeric MLE −101981.5; the β that maximises the closed-form likelihood is 0.50. -Optional exact-MLE refinement (argument `refine`, default 0): n fixed-point steps -wₜ = Qₜ^(β−1), α² = (β/T) Σ wₜ|eₜ|², ψ² = (β/T) Σ wₜ eₜ², then recompute ℓ with the full -density. Still closed-form arithmetic; the optimiser sees only β. Convergence is proven for -β ≤ 1 (Pascal et al., 2013); for β > 1 cap the iterations and fall back to step 2. +Exact MLE of ρ (fixed point wₜ = Qₜ^(β−1), α² = (β/T) Σ wₜ|eₜ|², ψ² = (β/T) Σ wₜ eₜ²) is NOT +implemented (DECIDED). Document in `?clm` and `?dcgnorm` that the scale/pseudo-scale estimates are +approximate ML for β ≠ 1: ρ is the moment estimator, the scale is exact ML given ρ. ## 3. Implementation phases @@ -55,22 +55,23 @@ Each phase ends with `R CMD build` + `R CMD check --no-manual` clean and its tes Commit per phase. ### Phase 1 — distribution functions (`R/cgnorm.R`, new) -- `dcgnorm(q, mu=0, sigma2=1, varsigma2=0, shape=1, log=FALSE)`. Default shape 1 = normal. +- `dcgnorm(q, mu=0, scale=1, pseudoscale=0, shape=1, log=FALSE)`. Default shape 1 = normal. In greybox `dgnorm` shape 2 = normal (exponent on |x|, not on the quadratic form); document that shape here equals greybox shape / 2. -- `rcgnorm(n, mu, sigma2, varsigma2, shape)` via the Gamma radial law and a uniform angle. -- Internal helper `cgnormScale(sigma2, varsigma2, shape)` returning (α², ψ²), and its inverse. -- Validate: sigma2 > 0, |varsigma2| < sigma2, shape > 0. +- `rcgnorm(n, mu, scale, pseudoscale, shape)` via the Gamma radial law and a uniform angle. +- Exported helper or documented formula k(β) to convert scale to covariance. +- Validate: scale > 0, |pseudoscale| < scale, shape > 0. - No `p`/`q` functions in this release (no closed form; `qcnorm` semantics are unusual anyway). - Roxygen docs in the `cnormal` style; export; NAMESPACE via roxygen. ### Phase 2 — `clm()` estimation (`R/clm.R`) -- New argument `distribution = c("dcnorm", "dcgnorm")`, used only when `loss = "likelihood"`. - Default keeps current behaviour exactly. -- Shape: pass via `...` as `shape` (NULL = estimate, number = fix), mirroring greybox `alm()` - handling of `other`. Returned in `object$other$shape`. +- New argument `distribution = c("dcnorm", "dcgnorm")` (DECIDED), as in greybox `alm()`: the + distribution defines the likelihood used when `loss = "likelihood"`. Default keeps current + behaviour exactly. +- Shape: as `alm()` does for `dgnorm`: user may fix it via `...` (`shape=`), otherwise estimated. + Returned in `object$other$shape`. - `fitter()` (≈ lines 195–247): for `dcgnorm`, compute ρ̂, α̂², ψ̂² as in section 2 and return - `scale` as the 2×2 covariance (via k(β)), so downstream code keeps its meaning. + them as `scale`. Check every downstream use of `object$scale` (vcov refit passes it back in). - Parameter vector: append ONE real parameter log β at the end of `BReal` in `estimator()` (≈ lines 307–341) and strip it BEFORE the real/imaginary split in `fitter()` (lines 197–201), which currently assumes every entry is half of a complex coefficient. Update `maxeval`. @@ -111,12 +112,38 @@ Commit per phase. 7. Regression: every existing example in the docs gives unchanged results with the default `distribution`. -## 5. Open questions for the author - -1. Variance parameterisation (recommended) or scale parameterisation for `dcgnorm`? -2. `distribution` as a new argument (recommended) or extend `loss` values (e.g. `"likelihood-cgn"`)? -3. Parameter counting: `nparam` counts complex units (one complex coefficient = 1, Σ = 1.5). - AIC/AICc of `clm` therefore penalise half the real-parameter count. Keep the convention (and - document it) or switch to real counts? This affects comparing cARIMA with real VAR by IC. -4. Ship `refine` (exact MLE of ρ) in the first version or later? -5. Name: `dcgnorm` (parallel to greybox `dgnorm`) — confirm. +## 5. Information criteria: parameter count (found while planning, affects master too) + +Check (cAR(1) example from the monograph, T = 80): +- `logLik.clm` is the FULL bivariate log-likelihood, −T(log 2π + 1 + ½ log det Σ̂) (−389.60). +- `nparam.clm` returns `object$rank` = `nVariables + 3/2` = 3.5, i.e. real parameters / 2 + (2 intercept + 2 AR + 3 Σ = 7 real). Commented-out code in `nparam.clm` says + "Divide by two, because we calculate df per series". +- `AIC`/`BIC` (stats, via `logLik` df attr) and `AICc`/`BICc` (greybox `.default`, via `nparam`) + all use 3.5. AIC = −2ℓ + 7 = 786.21; with all 7 real parameters it is −2ℓ + 14 = 793.21. +- For comparison greybox's own `AICc.varest` counts ALL real parameters of a VAR in the penalty, + K·m + K(K+1)/2 (coefficients of all equations plus Σ), and uses the per-equation count m only in + the small-sample denominator (Bedrick & Tsai 1994): −2ℓ + 2T(Km + K(K+1)/2)/(T − m − K − 1). + (`vars::logLik.varest` reports df = 6 for this VAR(1), omitting Σ; do not compare to it blindly.) + +Conclusion: the per-series count (3.5) is right for `df.residual` (each series has T obs and +3.5 parameters per series), but with a joint bivariate log-likelihood the AIC/BIC penalty needs +the total real count (7). The current values under-penalise by half, so AIC-based order +selection favours bigger cARIMA models, and IC comparison with real-valued VAR favours cARIMA. + +Proposed fix (Phase 0, separate commit, NEWS entry as a behaviour change): +- Keep `nparam.clm` / `df.residual` per series (no change to t-based inference). +- `logLik.clm`: `df = 2*nparam(object)` so stats `AIC`/`BIC` use the total real count. +- New `AICc.clm`, `BICc.clm` in the `AICc.varest` style with K = 2, k = 2*nparam(object) and + m = nparam(object) − (Σ and shape share) in the denominator. The exact small-sample + correction for the complex (restricted) model is not derived anywhere; decide the form of m. +- dcgnorm: β adds 1 real parameter (0.5 in the per-series count). +- Re-run the monograph examples: selected orders may change (BJsales chose cARIMA(3,1,1)). + +## 6. Decisions (2026-09-25) + +1. Scale parameterisation (`scale`, `pseudoscale`). +2. `distribution` argument, as in `alm()`. +3. IC count: see section 5, fix proposed, awaiting approval of the AICc/BICc form. +4. No exact-MLE refinement; document the approximation. +5. Name `dcgnorm` / `rcgnorm`. From b2389a2c6b4ce781af49274d6c046d142e9e5be5 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 09:44:18 +0000 Subject: [PATCH 03/16] Count all real parameters in AIC/BIC of clm; add multivariate AICc/BICc logLik.clm now reports df = 2*nparam (both parts of the complex variable), so AIC and BIC use the total number of real parameters. New AICc.clm and BICc.clm use the Bedrick & Tsai (1994) correction with two series, as in legion. nparam and df.residual remain per series. Adds testthat tests. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NAMESPACE | 2 ++ NEWS | 9 +++++++++ R/clm.R | 42 +++++++++++++++++++++++++++++++++++++++- dev/PLAN-dcgnorm.md | 4 +++- tests/testthat.R | 4 ++++ tests/testthat/test-ic.R | 38 ++++++++++++++++++++++++++++++++++++ 6 files changed, 97 insertions(+), 2 deletions(-) create mode 100644 tests/testthat.R create mode 100644 tests/testthat/test-ic.R diff --git a/NAMESPACE b/NAMESPACE index 7d94243..09c19d5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,7 @@ # Generated by roxygen2: do not edit by hand +S3method(AICc,clm) +S3method(BICc,clm) S3method(actuals,clm) S3method(confint,clm) S3method(logLik,clm) diff --git a/NEWS b/NEWS index 4eaee1d..2fe5caa 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,12 @@ +complex v1.0.3 (development) +============== + +Changes: +* AIC() and BIC() of clm() now count all real parameters of the model (both parts of the complex variable), as done in multivariate models. nparam() and the degrees of freedom are still per series. Before, the penalty was half of what it should be. Information criteria values change, and so can the orders selected by them. +* AICc() and BICc() methods for clm() use the multivariate small-sample correction (Bedrick & Tsai, 1994) with two series, as in the legion package. The correction assumes unrestricted multivariate regression, so it is approximate for complex models. +* Added tests (testthat). + + complex v1.0.2 (Release data: 2026-02-03) ============== diff --git a/R/clm.R b/R/clm.R index 95f2a87..f5248d5 100644 --- a/R/clm.R +++ b/R/clm.R @@ -994,7 +994,47 @@ nparam.clm <- function(object, all=TRUE, ...){ #' @importFrom stats logLik #' @export logLik.clm <- function(object, ...){ - return(structure(object$logLik,nobs=nobs(object),df=nparam(object),class="logLik")); + # nparam() counts per series (complex units). The likelihood is the joint one of the + # real and imaginary parts, so AIC/BIC need the total number of real parameters. + return(structure(object$logLik,nobs=nobs(object),df=2*nparam(object),class="logLik")); +} + +# Multivariate small-sample corrections (Bedrick & Tsai, 1994), as in legion, +# with two series (real and imaginary parts). nParamPerSeries excludes the share of +# the covariance matrix (3 real parameters, i.e. 3/2 per series). The correction is +# derived for unrestricted multivariate regression, so it is approximate for clm. +#' @export +AICc.clm <- function(object, ...){ + llikelihood <- logLik(object); + llikelihood <- llikelihood[1:length(llikelihood)]; + nSeries <- 2; + nParamAll <- 2*nparam(object); + nParamPerSeries <- nparam(object) - (nSeries+1)/2; + obs <- nobs(object); + if(obs - (nParamPerSeries + nSeries + 1) <= 0){ + IC <- Inf; + } + else{ + IC <- -2*llikelihood + 2*(obs*nParamAll/(obs - (nParamPerSeries + nSeries + 1))); + } + return(IC); +} + +#' @export +BICc.clm <- function(object, ...){ + llikelihood <- logLik(object); + llikelihood <- llikelihood[1:length(llikelihood)]; + nSeries <- 2; + nParamAll <- 2*nparam(object); + nParamPerSeries <- nparam(object) - (nSeries+1)/2; + obs <- nobs(object); + if(obs - (nParamPerSeries + nSeries + 1) <= 0){ + IC <- Inf; + } + else{ + IC <- -2*llikelihood + log(obs)*(obs*nParamAll/(obs - (nParamPerSeries + nSeries + 1))); + } + return(IC); } #' @rdname clm diff --git a/dev/PLAN-dcgnorm.md b/dev/PLAN-dcgnorm.md index a0ea13b..327c55c 100644 --- a/dev/PLAN-dcgnorm.md +++ b/dev/PLAN-dcgnorm.md @@ -144,6 +144,8 @@ Proposed fix (Phase 0, separate commit, NEWS entry as a behaviour change): 1. Scale parameterisation (`scale`, `pseudoscale`). 2. `distribution` argument, as in `alm()`. -3. IC count: see section 5, fix proposed, awaiting approval of the AICc/BICc form. +3. IC count: approved and implemented as proposed in section 5 (Phase 0 done, per-series nparam kept, + legion-style AICc/BICc). BJsales check: top-3 orders by AICc and BIC unchanged before/after + (cARIMA(3,1,2) first; the monograph reported (3,1,1), which is due to package version, not this fix). 4. No exact-MLE refinement; document the approximation. 5. Name `dcgnorm` / `rcgnorm`. diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..faf206a --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(complex) + +test_check("complex") diff --git a/tests/testthat/test-ic.R b/tests/testthat/test-ic.R new file mode 100644 index 0000000..fd38df1 --- /dev/null +++ b/tests/testthat/test-ic.R @@ -0,0 +1,38 @@ +# Information criteria of clm use the total number of real parameters +# (both parts of the complex variable), see Bedrick & Tsai (1994). +set.seed(41) +obs <- 110 +e <- rcnorm(obs, 0, sigma2=25, varsigma2=16+9i) +y <- vector("complex", obs) +y[1] <- 100+50i + e[1] +for(i in 2:obs){ + y[i] <- (0.5+0.2i)*y[i-1] + 100+50i + e[i] +} +y <- y[-c(1:10)] +cAR1Model <- clm(y~1, orders=c(1,0,0), subset=c(1:80)) + +test_that("nparam counts per series, logLik df counts all real parameters", { + # Intercept + AR(1): 2 complex = 4 real, covariance matrix: 3 real + expect_equal(nparam(cAR1Model), 3.5) + expect_equal(attr(logLik(cAR1Model), "df"), 7) + expect_equal(cAR1Model$df.residual, 80 - 3.5) +}) + +test_that("AIC and BIC use the total number of real parameters", { + ll <- as.numeric(logLik(cAR1Model)) + expect_equal(AIC(cAR1Model), -2*ll + 2*7) + expect_equal(BIC(cAR1Model), -2*ll + log(80)*7) +}) + +test_that("AICc and BICc use the multivariate correction with two series", { + ll <- as.numeric(logLik(cAR1Model)) + # 2 complex coefficients per series, 2 series + denominator <- 80 - (2 + 2 + 1) + expect_equal(AICc(cAR1Model), -2*ll + 2*80*7/denominator) + expect_equal(BICc(cAR1Model), -2*ll + log(80)*80*7/denominator) +}) + +test_that("ICs are not produced for non-likelihood losses", { + cAR1OLS <- clm(y~1, orders=c(1,0,0), subset=c(1:80), loss="OLS") + expect_true(is.na(logLik(cAR1OLS))) +}) From 4fbdc96f9b8493f144785b135d7e3c8bcfb52eae Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 09:44:24 +0000 Subject: [PATCH 04/16] Plan: correct note on BJsales order difference Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- dev/PLAN-dcgnorm.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/dev/PLAN-dcgnorm.md b/dev/PLAN-dcgnorm.md index 327c55c..928a18c 100644 --- a/dev/PLAN-dcgnorm.md +++ b/dev/PLAN-dcgnorm.md @@ -146,6 +146,7 @@ Proposed fix (Phase 0, separate commit, NEWS entry as a behaviour change): 2. `distribution` argument, as in `alm()`. 3. IC count: approved and implemented as proposed in section 5 (Phase 0 done, per-series nparam kept, legion-style AICc/BICc). BJsales check: top-3 orders by AICc and BIC unchanged before/after - (cARIMA(3,1,2) first; the monograph reported (3,1,1), which is due to package version, not this fix). + (cARIMA(3,1,2) first. The monograph reported (3,1,1); the old formula on the same fits also + gives (3,1,2), so the difference predates this fix. Cause not investigated.) 4. No exact-MLE refinement; document the approximation. 5. Name `dcgnorm` / `rcgnorm`. From ecb4120e49976c09e045d6dd2e8254a542b39353 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 09:57:14 +0000 Subject: [PATCH 05/16] Bump version to 1.0.3, regenerate NAMESPACE with roxygen2 8.1.0 Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- DESCRIPTION | 6 +-- NAMESPACE | 124 ++++++++++++++++++++++++++++------------------------ NEWS | 2 +- 3 files changed, 72 insertions(+), 60 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 67e1335..704baf9 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: complex Type: Package Title: Time Series Analysis and Forecasting Using Complex Variables -Version: 1.0.2 -Date: 2026-02-03 +Version: 1.0.3 +Date: 2026-09-25 Authors@R: person("Ivan", "Svetunkov", email = "ivan@svetunkov.com", role = c("aut", "cre"), comment="Lecturer at Centre for Marketing Analytics and Forecasting, Lancaster University, UK") URL: https://github.com/config-i1/complex @@ -27,6 +27,6 @@ Suggests: testthat, knitr, rmarkdown -RoxygenNote: 7.3.2 Encoding: UTF-8 Roxygen: list(old_usage = TRUE) +Config/roxygen2/version: 8.1.0 diff --git a/NAMESPACE b/NAMESPACE index 09c19d5..ea5553f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -39,62 +39,74 @@ export(qcnorm) export(rcnorm) export(vec2complex) import(legion) -importFrom(grDevices,devAskNewPage) -importFrom(grDevices,grey) -importFrom(grDevices,palette) -importFrom(graphics,abline) -importFrom(graphics,axis) -importFrom(graphics,box) -importFrom(graphics,layout) -importFrom(graphics,lines) -importFrom(graphics,mtext) -importFrom(graphics,par) -importFrom(graphics,points) -importFrom(graphics,text) -importFrom(greybox,AICc) -importFrom(greybox,BICc) -importFrom(greybox,actuals) -importFrom(greybox,measures) -importFrom(greybox,nparam) -importFrom(greybox,xregExpander) -importFrom(mvtnorm,dmvnorm) -importFrom(mvtnorm,pmvnorm) -importFrom(mvtnorm,qmvnorm) -importFrom(mvtnorm,rmvnorm) +importFrom(grDevices, + devAskNewPage, + grey, + palette +) +importFrom(graphics, + abline, + axis, + box, + layout, + lines, + mtext, + par, + points, + text +) +importFrom(greybox, + AICc, + BICc, + actuals, + measures, + nparam, + xregExpander +) +importFrom(mvtnorm, + dmvnorm, + pmvnorm, + qmvnorm, + rmvnorm +) importFrom(nloptr,nloptr) importFrom(pracma,hessian) -importFrom(stats,"contrasts<-") -importFrom(stats,AIC) -importFrom(stats,BIC) -importFrom(stats,as.formula) -importFrom(stats,cmdscale) -importFrom(stats,coef) -importFrom(stats,confint) -importFrom(stats,cor) -importFrom(stats,cov) -importFrom(stats,deltat) -importFrom(stats,dist) -importFrom(stats,fitted) -importFrom(stats,formula) -importFrom(stats,frequency) -importFrom(stats,logLik) -importFrom(stats,model.frame) -importFrom(stats,model.matrix) -importFrom(stats,nobs) -importFrom(stats,pacf) -importFrom(stats,qt) -importFrom(stats,resid) -importFrom(stats,residuals) -importFrom(stats,rnorm) -importFrom(stats,sd) -importFrom(stats,setNames) -importFrom(stats,sigma) -importFrom(stats,start) -importFrom(stats,terms) -importFrom(stats,time) -importFrom(stats,ts) -importFrom(stats,update.formula) -importFrom(stats,vcov) -importFrom(utils,packageVersion) -importFrom(utils,tail) +importFrom(stats, + "contrasts<-", + AIC, + BIC, + as.formula, + cmdscale, + coef, + confint, + cor, + cov, + deltat, + dist, + fitted, + formula, + frequency, + logLik, + model.frame, + model.matrix, + nobs, + pacf, + qt, + resid, + residuals, + rnorm, + sd, + setNames, + sigma, + start, + terms, + time, + ts, + update.formula, + vcov +) +importFrom(utils, + packageVersion, + tail +) useDynLib(complex) diff --git a/NEWS b/NEWS index 2fe5caa..6ff0c3c 100644 --- a/NEWS +++ b/NEWS @@ -1,4 +1,4 @@ -complex v1.0.3 (development) +complex v1.0.3 (Release data: 2026-09-25) ============== Changes: From 5c1a21ec61d4629520ff38f440f40dadb75f9ae3 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 10:02:28 +0000 Subject: [PATCH 06/16] Update repository links to openforecast-org/complex Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- DESCRIPTION | 7 ++++--- README.md | 6 +++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 704baf9..42de076 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -4,9 +4,10 @@ Title: Time Series Analysis and Forecasting Using Complex Variables Version: 1.0.3 Date: 2026-09-25 Authors@R: person("Ivan", "Svetunkov", email = "ivan@svetunkov.com", role = c("aut", "cre"), - comment="Lecturer at Centre for Marketing Analytics and Forecasting, Lancaster University, UK") -URL: https://github.com/config-i1/complex -BugReports: https://github.com/config-i1/complex/issues + comment = c("Director of the OpenForecast consultancy company, UK", + ORCID = "0000-0001-7826-0281")) +URL: https://github.com/openforecast-org/complex +BugReports: https://github.com/openforecast-org/complex/issues Language: en-GB Description: Implements the instruments for complex-valued modelling, including time series analysis and forecasting. This is based on the monograph diff --git a/README.md b/README.md index ccc18a2..895f85f 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,11 @@ # complex package for R [![CRAN_Status_Badge](https://www.r-pkg.org/badges/version/complex)](https://cran.r-project.org/package=complex) [![Downloads](https://cranlogs.r-pkg.org/badges/complex)](https://cran.r-project.org/package=complex) -[![R-CMD-check](https://github.com/config-i1/complex/actions/workflows/test.yml/badge.svg)](https://github.com/config-i1/complex/actions/workflows/test.yml) +[![R-CMD-check](https://github.com/openforecast-org/complex/actions/workflows/test.yml/badge.svg)](https://github.com/openforecast-org/complex/actions/workflows/test.yml) Time series analysis and forecasting using complex variables -![hex-sticker of the complex package for R](https://github.com/config-i1/complex/blob/master/man/figures/complex-web.png?raw=true) +![hex-sticker of the complex package for R](https://github.com/openforecast-org/complex/blob/master/man/figures/complex-web.png?raw=true) The package includes basic instruments for correlation and regression analysis of complex-valued variables. The package supports the monograph by Svetunkov & Svetunkov "Complex-valued Econometrics with Examples in R", which is to be published by Springer in 2024. @@ -19,5 +19,5 @@ install.packages("complex") For installation from github use the `remotes` package: ```r if (!require("remotes")){install.packages("remotes")} -remotes::install_github("config-i1/complex") +remotes::install_github("openforecast-org/complex") ``` From 3158b542a2676f6df801a93cc3e8c3537402fa61 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 14:10:12 +0000 Subject: [PATCH 07/16] Add complex generalised normal distribution (dcgnorm) and use it in clm() dcgnorm()/rcgnorm() in scale parameterisation. clm() gets a distribution argument; for dcgnorm the shape is estimated (or fixed via shape=) and the scale/pseudo-scale are concentrated out in closed form. nparam counts the shape, AICc/BICc exclude it from the per-series regressors, summary prints it. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NAMESPACE | 4 ++ NEWS | 2 + R/cgnorm.R | 111 +++++++++++++++++++++++++++++++++++ R/clm.R | 86 +++++++++++++++++++++++++-- dev/PLAN-dcgnorm.md | 10 ++++ man/cgnormal.Rd | 81 +++++++++++++++++++++++++ man/clm.Rd | 16 ++++- tests/testthat/test-cgnorm.R | 82 ++++++++++++++++++++++++++ 8 files changed, 384 insertions(+), 8 deletions(-) create mode 100644 R/cgnorm.R create mode 100644 man/cgnormal.Rd create mode 100644 tests/testthat/test-cgnorm.R diff --git a/NAMESPACE b/NAMESPACE index ea5553f..896a87a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -31,11 +31,13 @@ export(cpacf) export(cplot) export(cscale) export(cvar) +export(dcgnorm) export(dcnorm) export(invert) export(mat2complex) export(pcnorm) export(qcnorm) +export(rcgnorm) export(rcnorm) export(vec2complex) import(legion) @@ -94,7 +96,9 @@ importFrom(stats, qt, resid, residuals, + rgamma, rnorm, + runif, sd, setNames, sigma, diff --git a/NEWS b/NEWS index 6ff0c3c..0f70021 100644 --- a/NEWS +++ b/NEWS @@ -2,6 +2,8 @@ complex v1.0.3 (Release data: 2026-09-25) ============== Changes: +* dcgnorm() and rcgnorm() - density and random generation for the complex generalised normal distribution (complex generalised Gaussian), with scale, pseudo-scale and shape. +* clm() has a new argument distribution. distribution="dcgnorm" uses the complex generalised normal likelihood, with the shape estimated (or provided via shape=...). The scale and the pseudo-scale are concentrated out: the circularity coefficient comes from the moments of the residuals and the scale is the ML estimate given it (approximate ML when shape is not 1). Standard errors are for the parameters given the shape. * AIC() and BIC() of clm() now count all real parameters of the model (both parts of the complex variable), as done in multivariate models. nparam() and the degrees of freedom are still per series. Before, the penalty was half of what it should be. Information criteria values change, and so can the orders selected by them. * AICc() and BICc() methods for clm() use the multivariate small-sample correction (Bedrick & Tsai, 1994) with two series, as in the legion package. The correction assumes unrestricted multivariate regression, so it is approximate for complex models. * Added tests (testthat). diff --git a/R/cgnorm.R b/R/cgnorm.R new file mode 100644 index 0000000..e7a66bb --- /dev/null +++ b/R/cgnorm.R @@ -0,0 +1,111 @@ +#' Complex Generalised Normal Distribution +#' +#' Density and random generation for the complex generalised normal distribution +#' (also known as the complex generalised Gaussian or the bivariate exponential power +#' distribution) with location \code{mu}, scale \code{scale}, pseudo-scale +#' \code{pseudoscale} and shape \code{shape}. +#' +#' For the error term \eqn{e = z - \mu}, the density is +#' \deqn{f(z) = \frac{2\beta}{\pi \Gamma(1/\beta) 2^{1/\beta} \sqrt{\alpha^4 - |\psi^2|^2}} +#' \exp\left(-\frac{1}{2} Q^\beta\right),} +#' where \eqn{Q = 2(\alpha^2 |e|^2 - \Re(\overline{\psi^2} e^2)) / (\alpha^4 - |\psi^2|^2)}, +#' \eqn{\alpha^2} is the scale, \eqn{\psi^2} is the pseudo-scale and \eqn{\beta} is the shape. +#' +#' With \code{shape=1} this is the complex normal distribution with variance \code{scale} and +#' pseudo-variance \code{pseudoscale} (see \link[complex]{dcnorm}). Values of shape below one +#' give heavier tails (\code{shape=0.5} is Laplace-like), values above one give lighter tails +#' (the uniform distribution on an ellipse in the limit). Note that the shape here applies to +#' the quadratic form, so \code{shape=1} corresponds to \code{shape=2} in +#' \link[greybox]{dgnorm}. +#' +#' The variance and the pseudo-variance of the distribution equal \code{scale} and +#' \code{pseudoscale} multiplied by \eqn{2^{1/\beta} \Gamma(2/\beta) / (2 \Gamma(1/\beta))}. +#' +#' Random numbers are generated using the fact that \eqn{Q^\beta} follows the Gamma distribution +#' with shape \eqn{1/\beta} and scale 2, and the angle is uniform in the standardised space. +#' +#' @template author +#' @keywords distribution +#' +#' @param q vector of quantiles (complex). +#' @param n number of observations. Should be a single number. +#' @param mu vector of location parameters (complex). +#' @param scale scale parameter (real, positive). +#' @param pseudoscale pseudo-scale parameter (complex, with the modulus smaller than scale). +#' @param shape shape parameter (real, positive). +#' @param log if \code{TRUE}, then the logarithm of the density is returned. +#' +#' @return \code{dcgnorm} returns the density, \code{rcgnorm} returns a vector of +#' complex random values. +#' +#' @references \itemize{ +#' \item Novey, M., Adali, T., Roy, A. (2010). A complex generalized Gaussian distribution - +#' characterization, generation, and estimation. IEEE Transactions on Signal Processing, +#' 58(3), 1427-1433. +#' \item Gomez, E., Gomez-Villegas, M. A., Marin, J. M. (1998). A multivariate generalization +#' of the power exponential family of distributions. Communications in Statistics - Theory +#' and Methods, 27(3), 589-600. +#' } +#' +#' @seealso \code{\link[complex]{dcnorm}} +#' +#' @examples +#' x <- rcgnorm(1000, 1+1i, scale=2, pseudoscale=0.6+0.8i, shape=0.5) +#' plot(x) +#' dcgnorm(x[1:5], 1+1i, scale=2, pseudoscale=0.6+0.8i, shape=0.5) +#' +#' @rdname cgnormal +#' @aliases cgnormal dcgnorm +#' @export dcgnorm +dcgnorm <- function(q, mu=0, scale=1, pseudoscale=0, shape=1, log=FALSE){ + cgnormCheck(scale, pseudoscale, shape); + errors <- as.vector(q - mu); + determinant <- scale^2 - Mod(pseudoscale)^2; + quadraticForm <- 2*(scale*Mod(errors)^2 - Re(Conj(pseudoscale)*errors^2)) / determinant; + cgnormReturn <- (log(2*shape) - log(pi) - lgamma(1/shape) - log(2)/shape - + 0.5*log(determinant) - 0.5*quadraticForm^shape); + if(!log){ + cgnormReturn <- exp(cgnormReturn); + } + return(cgnormReturn); +} + +#' @rdname cgnormal +#' @export rcgnorm +#' @aliases rcgnorm +#' @importFrom stats rgamma runif +rcgnorm <- function(n=1, mu=0, scale=1, pseudoscale=0, shape=1){ + cgnormCheck(scale, pseudoscale, shape); + # Covariance-type matrix of the real and imaginary parts (as in dcnorm) + Sigma <- matrix(c(scale+Re(pseudoscale), Im(pseudoscale), + Im(pseudoscale), scale-Re(pseudoscale))/2, 2, 2); + radius <- sqrt(rgamma(n, shape=1/shape, scale=2)^(1/shape)); + angle <- runif(n, 0, 2*pi); + cgnormReturn <- t(chol(Sigma)) %*% rbind(radius*cos(angle), radius*sin(angle)); + return(mu + complex(real=cgnormReturn[1,], imaginary=cgnormReturn[2,])); +} + +cgnormCheck <- function(scale, pseudoscale, shape){ + if(any(scale<=0)){ + stop("The scale should be positive.", call.=FALSE); + } + if(any(Mod(pseudoscale)>=scale)){ + stop("The modulus of the pseudoscale should be smaller than the scale.", call.=FALSE); + } + if(any(shape<=0)){ + stop("The shape should be positive.", call.=FALSE); + } +} + +# Concentrated likelihood of the complex generalised normal for given errors and shape. +# The circularity coefficient rho = pseudoscale / scale is estimated from moments +# (exact ML only when shape=1); the scale is the exact ML estimate given rho. +cgnormConcentrated <- function(errors, shape){ + obs <- length(errors); + rho <- sum(errors^2) / sum(Mod(errors)^2); + quadraticForm <- 2*(Mod(errors)^2 - Re(Conj(rho)*errors^2)) / (1-Mod(rho)^2); + scale <- (shape/(2*obs) * sum(quadraticForm^shape))^(1/shape); + logLik <- obs*(log(2*shape) - log(pi) - lgamma(1/shape) - log(2)/shape - + log(scale) - 0.5*log(1-Mod(rho)^2) - 1/shape); + return(list(scale=scale, pseudoscale=rho*scale, logLik=logLik)); +} diff --git a/R/clm.R b/R/clm.R index f5248d5..a5b2f89 100644 --- a/R/clm.R +++ b/R/clm.R @@ -41,6 +41,16 @@ #' \code{lossFunction <- function(actual, fitted, B, xreg) return(mean(abs(actual-fitted)))} #' \code{loss=lossFunction} #' +#' @param distribution the distribution used in the likelihood, when \code{loss="likelihood"}: +#' \itemize{ +#' \item \code{"dcnorm"} - complex normal distribution (see \link[complex]{dcnorm}); +#' \item \code{"dcgnorm"} - complex generalised normal distribution +#' (see \link[complex]{dcgnorm}). Its shape parameter is estimated, unless it is provided via +#' \code{shape} in the ellipsis. The scale and the pseudo-scale are not estimated by the optimiser: +#' the circularity coefficient (pseudo-scale divided by scale) is taken from the moments of +#' the residuals and the scale is the maximum likelihood estimate given that coefficient. This is +#' exact maximum likelihood only for \code{shape=1}, and an approximation otherwise. +#' } #' @param orders vector of orders of complex ARIMA(p,d,q). #' @param scaling NOT YET IMPLEMENTED!!! Defines what type of scaling to do for the variables. #' See \link[complex]{cscale} for the explanation of the options. @@ -55,6 +65,8 @@ #' \item \code{FI=TRUE} will make the function also produce Fisher Information #' matrix, which then can be used to calculated variances of smoothing parameters #' and initial states of the model. This is used in the \link[stats]{vcov} method; +#' \item \code{shape} - the value of the shape parameter of \code{distribution="dcgnorm"}. If +#' provided, the shape is not estimated; #' } #' #' You can also pass parameters to the optimiser: @@ -136,6 +148,7 @@ #' @export clm clm <- function(formula, data, subset, na.action, loss=c("likelihood","OLS","CLS","MSE","MAE","HAM"), + distribution=c("dcnorm","dcgnorm"), orders=c(0,0,0), scaling=c("normalisation","standardisation","max","none"), parameters=NULL, fast=FALSE, ...){ # Start measuring the time of calculations @@ -157,6 +170,10 @@ clm <- function(formula, data, subset, na.action, } scaling <- match.arg(scaling); + distribution <- match.arg(distribution); + if(distribution=="dcgnorm" && loss!="likelihood"){ + stop("distribution=\"dcgnorm\" can only be used with loss=\"likelihood\".", call.=FALSE); + } #### Functions used in the estimation #### ifelseFast <- function(condition, yes, no){ @@ -257,9 +274,18 @@ clm <- function(formula, data, subset, na.action, } CF <- function(B, loss, y, matrixXreg){ + # If the shape is estimated, its logarithm is the last element of the real-valued vector + shapeValue <- shape; + if(shapeEstimate && !is.complex(B) && length(B) %% 2 == 1){ + shapeValue <- exp(B[length(B)]); + B <- B[-length(B)]; + } fitterReturn <- fitter(B, y, matrixXreg); - if(loss=="likelihood"){ + if(loss=="likelihood" && distribution=="dcgnorm"){ + CFValue <- -cgnormConcentrated(y - fitterReturn$mu, shapeValue)$logLik; + } + else if(loss=="likelihood"){ # # Concentrated logLik CFValue <- obsInsample*(log(2*pi) + 1 + 0.5*log(det(fitterReturn$scale))); @@ -321,6 +347,11 @@ clm <- function(formula, data, subset, na.action, # Retransform the vector of parameters into a real-valued one BReal <- c(Re(B),Im(B)); nVariables <- length(BReal); + # The logarithm of the shape goes to the end of the vector + if(shapeEstimate){ + BReal <- c(BReal, log(shape)); + maxeval <- maxeval + 40; + } # Although this is not needed in case of distribution="dnorm", we do that in a way, for the code consistency purposes res <- nloptr(BReal, CF, @@ -329,6 +360,10 @@ clm <- function(formula, data, subset, na.action, # lb=BLower, ub=BUpper, loss=loss, y=y, matrixXreg=matrixXreg); BReal[] <- res$solution; + shapeValue <- shape; + if(shapeEstimate){ + shapeValue <- exp(BReal[nVariables+1]); + } B[] <- complex(real=BReal[1:(nVariables/2)],imaginary=BReal[(nVariables/2+1):nVariables]); nVariables <- length(B); CFValue <- res$objective; @@ -337,13 +372,32 @@ clm <- function(formula, data, subset, na.action, print(res); } - return(list(B=B, CFValue=CFValue)); + return(list(B=B, CFValue=CFValue, shape=shapeValue)); } #### Define the rest of parameters #### ellipsis <- list(...); # ellipsis <- match.call(expand.dots = FALSE)$`...`; + # Shape of the complex generalised normal: estimated unless provided + shape <- 1; + shapeEstimate <- FALSE; + if(distribution=="dcgnorm"){ + if(is.null(ellipsis$shape)){ + if(!is.null(parameters)){ + stop("When parameters are provided for distribution=\"dcgnorm\", ", + "the shape needs to be provided as well.", call.=FALSE); + } + shapeEstimate <- TRUE; + } + else{ + shape <- ellipsis$shape; + if(shape<=0){ + stop("The shape should be positive.", call.=FALSE); + } + } + } + # Fisher Information if(is.null(ellipsis$FI)){ FI <- FALSE; @@ -861,6 +915,7 @@ clm <- function(formula, data, subset, na.action, res <- estimator(B, print_level); B <- res$B; CFValue <- res$CFValue; + shape <- res$shape; } } # If the parameters are provided @@ -895,6 +950,12 @@ clm <- function(formula, data, subset, na.action, mu[] <- fitterReturn$mu; scale <- fitterReturn$scale; matrixXreg[] <- fitterReturn$matrixXreg; + if(distribution=="dcgnorm"){ + cgnormValues <- cgnormConcentrated(y - mu, shape); + scale <- c(scale=cgnormValues$scale, pseudoscale=cgnormValues$pseudoscale); + ellipsis$shape <- shape; + ellipsis$shapeEstimated <- shapeEstimate; + } #### Produce Fisher Information #### if(FI){ @@ -922,7 +983,8 @@ clm <- function(formula, data, subset, na.action, ### Error term in the transformed scale errors[] <- extractorResiduals(mu, yFitted); - nParam <- nVariables + (loss=="likelihood")*3/2; + # Parameters are counted per series: 3/2 for the covariance matrix, 1/2 for the shape + nParam <- nVariables + (loss=="likelihood")*3/2 + shapeEstimate/2; if(interceptIsNeeded){ # This shit is needed, because R has habit of converting everything into vectors... @@ -959,7 +1021,8 @@ clm <- function(formula, data, subset, na.action, #### Return the model #### finalModel <- structure(list(coefficients=parameters, FI=FI, fitted=yFitted, residuals=as.vector(errors), mu=mu, scale=scale, logLik=logLik, model=modelName, - loss=loss, lossFunction=lossFunction, lossValue=CFValue, + loss=loss, distribution=distribution, + lossFunction=lossFunction, lossValue=CFValue, df.residual=obsInsample-nParam, df=nParam, call=cl, rank=nParam, data=dataWork, terms=dataTerms, subset=subset, other=ellipsis, B=B, @@ -1009,7 +1072,7 @@ AICc.clm <- function(object, ...){ llikelihood <- llikelihood[1:length(llikelihood)]; nSeries <- 2; nParamAll <- 2*nparam(object); - nParamPerSeries <- nparam(object) - (nSeries+1)/2; + nParamPerSeries <- nparam(object) - (nSeries+1)/2 - isTRUE(object$other$shapeEstimated)/2; obs <- nobs(object); if(obs - (nParamPerSeries + nSeries + 1) <= 0){ IC <- Inf; @@ -1026,7 +1089,7 @@ BICc.clm <- function(object, ...){ llikelihood <- llikelihood[1:length(llikelihood)]; nSeries <- 2; nParamAll <- 2*nparam(object); - nParamPerSeries <- nparam(object) - (nSeries+1)/2; + nParamPerSeries <- nparam(object) - (nSeries+1)/2 - isTRUE(object$other$shapeEstimated)/2; obs <- nobs(object); if(obs - (nParamPerSeries + nSeries + 1) <= 0){ IC <- Inf; @@ -1213,6 +1276,11 @@ vcov.clm <- function(object, type=NULL, ...){ newCall$orders <- object$other$orders; newCall$parameters <- c(Re(coef(object)),Im(coef(object))); newCall$scale <- object$scale; + # The shape of dcgnorm is fixed at its estimate: the Hessian is for B given the shape + if(!is.null(object$distribution) && object$distribution=="dcgnorm"){ + newCall$distribution <- "dcgnorm"; + newCall$shape <- object$other$shape; + } newCall$fast <- TRUE; newCall$FI <- TRUE; # Include bloody ellipsis @@ -1334,6 +1402,7 @@ summary.clm <- function(object, level=0.95, ...){ ourReturn$ICs <- ICs; } ourReturn$loss <- object$loss; + ourReturn$distribution <- object$distribution; ourReturn$model <- object$model; ourReturn$other <- object$other; ourReturn$responseName <- formula(object)[[2]]; @@ -1365,6 +1434,11 @@ print.summary.clm <- function(x, ...){ cat(x$model, "estimated via clm()\n"); cat(paste0("Response variable: ", paste0(x$responseName,collapse=""))); cat(paste0("\nLoss function used in estimation: ",x$loss)); + if(!is.null(x$distribution) && x$distribution=="dcgnorm"){ + cat(paste0("\nDistribution: complex generalised normal, shape=", + round(x$other$shape,digits), + ifelse(isTRUE(x$other$shapeEstimated), " (estimated)", " (provided)"))); + } cat("\nCoefficients:\n"); stars <- setNames(vector("character",length(x$significance)), diff --git a/dev/PLAN-dcgnorm.md b/dev/PLAN-dcgnorm.md index 928a18c..a1d3db8 100644 --- a/dev/PLAN-dcgnorm.md +++ b/dev/PLAN-dcgnorm.md @@ -49,6 +49,16 @@ Exact MLE of ρ (fixed point wₜ = Qₜ^(β−1), α² = (β/T) Σ wₜ|eₜ|² implemented (DECIDED). Document in `?clm` and `?dcgnorm` that the scale/pseudo-scale estimates are approximate ML for β ≠ 1: ρ is the moment estimator, the scale is exact ML given ρ. +## Status (2026-09-25) + +Phase 0 done (IC fix). Phases 1 and 2 done: `R/cgnorm.R` (dcgnorm, rcgnorm, internal +cgnormConcentrated), clm(distribution=), shape estimated as log-shape appended to the optimiser +vector, nparam +0.5, AICc/BICc m excludes the shape, vcov refit fixes the shape (Hessian for B +given shape, Phase 3 still to do properly), summary prints the shape. Tests in +tests/testthat/test-cgnorm.R. Remaining: Phase 3 (joint Hessian), Phase 4 (simulated intervals), +Phase 5 (diagnostic plot), and clm's poor default starting values (found in the wind experiment: +nested bigger models can end with lower logLik). + ## 3. Implementation phases Each phase ends with `R CMD build` + `R CMD check --no-manual` clean and its tests passing. diff --git a/man/cgnormal.Rd b/man/cgnormal.Rd new file mode 100644 index 0000000..84d8650 --- /dev/null +++ b/man/cgnormal.Rd @@ -0,0 +1,81 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cgnorm.R +\name{dcgnorm} +\alias{dcgnorm} +\alias{cgnormal} +\alias{rcgnorm} +\title{Complex Generalised Normal Distribution} +\usage{ +dcgnorm(q, mu = 0, scale = 1, pseudoscale = 0, shape = 1, + log = FALSE) + +rcgnorm(n = 1, mu = 0, scale = 1, pseudoscale = 0, shape = 1) +} +\arguments{ +\item{q}{vector of quantiles (complex).} + +\item{mu}{vector of location parameters (complex).} + +\item{scale}{scale parameter (real, positive).} + +\item{pseudoscale}{pseudo-scale parameter (complex, with the modulus smaller than scale).} + +\item{shape}{shape parameter (real, positive).} + +\item{log}{if \code{TRUE}, then the logarithm of the density is returned.} + +\item{n}{number of observations. Should be a single number.} +} +\value{ +\code{dcgnorm} returns the density, \code{rcgnorm} returns a vector of +complex random values. +} +\description{ +Density and random generation for the complex generalised normal distribution +(also known as the complex generalised Gaussian or the bivariate exponential power +distribution) with location \code{mu}, scale \code{scale}, pseudo-scale +\code{pseudoscale} and shape \code{shape}. +} +\details{ +For the error term \eqn{e = z - \mu}, the density is +\deqn{f(z) = \frac{2\beta}{\pi \Gamma(1/\beta) 2^{1/\beta} \sqrt{\alpha^4 - |\psi^2|^2}} +\exp\left(-\frac{1}{2} Q^\beta\right),} +where \eqn{Q = 2(\alpha^2 |e|^2 - \Re(\overline{\psi^2} e^2)) / (\alpha^4 - |\psi^2|^2)}, +\eqn{\alpha^2} is the scale, \eqn{\psi^2} is the pseudo-scale and \eqn{\beta} is the shape. + +With \code{shape=1} this is the complex normal distribution with variance \code{scale} and +pseudo-variance \code{pseudoscale} (see \link[complex]{dcnorm}). Values of shape below one +give heavier tails (\code{shape=0.5} is Laplace-like), values above one give lighter tails +(the uniform distribution on an ellipse in the limit). Note that the shape here applies to +the quadratic form, so \code{shape=1} corresponds to \code{shape=2} in +\link[greybox]{dgnorm}. + +The variance and the pseudo-variance of the distribution equal \code{scale} and +\code{pseudoscale} multiplied by \eqn{2^{1/\beta} \Gamma(2/\beta) / (2 \Gamma(1/\beta))}. + +Random numbers are generated using the fact that \eqn{Q^\beta} follows the Gamma distribution +with shape \eqn{1/\beta} and scale 2, and the angle is uniform in the standardised space. +} +\examples{ +x <- rcgnorm(1000, 1+1i, scale=2, pseudoscale=0.6+0.8i, shape=0.5) +plot(x) +dcgnorm(x[1:5], 1+1i, scale=2, pseudoscale=0.6+0.8i, shape=0.5) + +} +\references{ +\itemize{ +\item Novey, M., Adali, T., Roy, A. (2010). A complex generalized Gaussian distribution - +characterization, generation, and estimation. IEEE Transactions on Signal Processing, +58(3), 1427-1433. +\item Gomez, E., Gomez-Villegas, M. A., Marin, J. M. (1998). A multivariate generalization +of the power exponential family of distributions. Communications in Statistics - Theory +and Methods, 27(3), 589-600. +} +} +\seealso{ +\code{\link[complex]{dcnorm}} +} +\author{ +Ivan Svetunkov, \email{ivan@svetunkov.com} +} +\keyword{distribution} diff --git a/man/clm.Rd b/man/clm.Rd index 207e0fe..7354580 100644 --- a/man/clm.Rd +++ b/man/clm.Rd @@ -8,8 +8,9 @@ \title{Complex Linear Model} \usage{ clm(formula, data, subset, na.action, loss = c("likelihood", "OLS", "CLS", - "MSE", "MAE", "HAM"), orders = c(0, 0, 0), scaling = c("normalisation", - "standardisation", "max", "none"), parameters = NULL, fast = FALSE, ...) + "MSE", "MAE", "HAM"), distribution = c("dcnorm", "dcgnorm"), + orders = c(0, 0, 0), scaling = c("normalisation", "standardisation", + "max", "none"), parameters = NULL, fast = FALSE, ...) \method{sigma}{clm}(object, type = NULL, ...) @@ -54,6 +55,17 @@ example: \code{lossFunction <- function(actual, fitted, B, xreg) return(mean(abs(actual-fitted)))} \code{loss=lossFunction}} +\item{distribution}{the distribution used in the likelihood, when \code{loss="likelihood"}: +\itemize{ +\item \code{"dcnorm"} - complex normal distribution (see \link[complex]{dcnorm}); +\item \code{"dcgnorm"} - complex generalised normal distribution +(see \link[complex]{dcgnorm}). Its shape parameter is estimated, unless it is provided via +\code{shape} in the ellipsis. The scale and the pseudo-scale are not estimated by the optimiser: +the circularity coefficient (pseudo-scale divided by scale) is taken from the moments of +the residuals and the scale is the maximum likelihood estimate given that coefficient. This is +exact maximum likelihood only for \code{shape=1}, and an approximation otherwise. +}} + \item{orders}{vector of orders of complex ARIMA(p,d,q).} \item{scaling}{NOT YET IMPLEMENTED!!! Defines what type of scaling to do for the variables. diff --git a/tests/testthat/test-cgnorm.R b/tests/testthat/test-cgnorm.R new file mode 100644 index 0000000..f311015 --- /dev/null +++ b/tests/testthat/test-cgnorm.R @@ -0,0 +1,82 @@ +# Complex generalised normal distribution and its use in clm() + +test_that("dcgnorm with shape=1 equals dcnorm", { + set.seed(41) + z <- rcnorm(20, 1+1i, sigma2=2, varsigma2=0.6+0.8i) + expect_equal(dcgnorm(z, 1+1i, scale=2, pseudoscale=0.6+0.8i, shape=1, log=TRUE), + dcnorm(z, 1+1i, sigma2=2, varsigma2=0.6+0.8i, log=TRUE), + tolerance=1e-12) +}) + +test_that("dcgnorm integrates to one", { + grid <- seq(-40, 40, length.out=801) + step <- diff(grid)[1] + z <- outer(grid, 1i*grid, "+") + for(shape in c(0.5, 1, 2)){ + expect_equal(sum(dcgnorm(z, 0, scale=2, pseudoscale=0.6+0.8i, shape=shape))*step^2, + 1, tolerance=1e-3) + } +}) + +test_that("rcgnorm has the right covariance and radial distribution", { + set.seed(41) + shape <- 0.5 + scale <- 2 + pseudoscale <- 0.6+0.8i + z <- rcgnorm(50000, 0, scale, pseudoscale, shape) + # Variance and pseudo-variance are the scale and pseudo-scale times this factor + multiplier <- 2^(1/shape)*gamma(2/shape)/(2*gamma(1/shape)) + expect_equal(mean(Mod(z)^2), multiplier*scale, tolerance=0.05) + expect_equal(mean(z^2), multiplier*pseudoscale, tolerance=0.05) + # Q^shape follows Gamma(1/shape, scale=2) + quadraticForm <- 2*(scale*Mod(z)^2 - Re(Conj(pseudoscale)*z^2))/(scale^2-Mod(pseudoscale)^2) + expect_gt(suppressWarnings(ks.test(quadraticForm^shape, "pgamma", shape=1/shape, scale=2)$p.value), + 0.001) +}) + +test_that("the concentrated likelihood matches the density at the estimated scales", { + set.seed(41) + e <- rcgnorm(2000, 0, 2, 0.6+0.8i, 0.5) + values <- complex:::cgnormConcentrated(e, 0.5) + expect_equal(values$logLik, + sum(dcgnorm(e, 0, values$scale, values$pseudoscale, 0.5, log=TRUE))) +}) + +# cAR(1) with heavy-tailed errors +set.seed(41) +obs <- 2100 +e <- rcgnorm(obs, 0, scale=4, pseudoscale=1+1i, shape=0.5) +y <- vector("complex", obs) +y[1] <- 10+5i + e[1] +for(i in 2:obs){ + y[i] <- (0.5+0.2i)*y[i-1] + 10+5i + e[i] +} +y <- y[-c(1:100)] +data <- data.frame(y=y) + +test_that("clm with dcgnorm and shape=1 reproduces dcnorm", { + modelNormal <- clm(y~1, data, orders=c(1,0,0)) + modelShapeOne <- clm(y~1, data, orders=c(1,0,0), distribution="dcgnorm", shape=1) + expect_equal(coef(modelShapeOne), coef(modelNormal), tolerance=1e-4) + expect_equal(as.numeric(logLik(modelShapeOne)), as.numeric(logLik(modelNormal)), + tolerance=1e-6) + expect_equal(nparam(modelShapeOne), nparam(modelNormal)) +}) + +test_that("clm with dcgnorm estimates the shape and beats dcnorm on heavy tails", { + modelNormal <- clm(y~1, data, orders=c(1,0,0)) + modelGeneral <- clm(y~1, data, orders=c(1,0,0), distribution="dcgnorm") + expect_true(modelGeneral$other$shapeEstimated) + expect_gt(modelGeneral$other$shape, 0.4) + expect_lt(modelGeneral$other$shape, 0.6) + expect_equal(nparam(modelGeneral), nparam(modelNormal) + 0.5) + expect_lt(AICc(modelGeneral), AICc(modelNormal)) + expect_equal(unname(coef(modelGeneral)[2]), 0.5+0.2i, tolerance=0.05) + # summary and forecasts work + expect_output(print(summary(modelGeneral)), "shape") + expect_length(predict(modelGeneral, newdata=matrix(NA, 5, 1))$mean, 5) +}) + +test_that("dcgnorm is refused with non-likelihood losses", { + expect_error(clm(y~1, data, loss="OLS", distribution="dcgnorm")) +}) From d4f8d7cb3b2885a9509564a1d7c356607e3b90d8 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 14:21:01 +0000 Subject: [PATCH 08/16] clm: accept a starting shape for dcgnorm as the last element of B; fix ... docs The @param ... of sigma.clm (same Rd page) overrode clm's documentation of the ellipsis (optimiser options, FI, shape). Removed it so that clm.Rd shows the full list. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- R/clm.R | 9 ++++++++- man/clm.Rd | 34 +++++++++++++++++++++++++++++++++- 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/R/clm.R b/R/clm.R index a5b2f89..1cec842 100644 --- a/R/clm.R +++ b/R/clm.R @@ -74,6 +74,8 @@ #' \item \code{B} - the vector of starting values of parameters for the optimiser, #' should correspond to the explanatory variables. If formula for scale was provided, #' the parameters for that part should follow the parameters for location; +#' for \code{distribution="dcgnorm"} with the estimated shape, the starting value of the shape +#' can be added as the last element of \code{B}; #' \item \code{algorithm} - the algorithm to use in optimisation #' (\code{"NLOPT_LN_SBPLX"} by default); #' \item \code{maxeval} - maximum number of evaluations to carry out. Default is 40 per @@ -912,6 +914,12 @@ clm <- function(formula, data, subset, na.action, B <- c(B, rep(0.1*(1+1i),maOrder)); } } + # For dcgnorm with the estimated shape, B can have one extra element in the end: + # the starting value of the shape (its real part) + if(shapeEstimate && length(B)==length(parametersNames)+1){ + shape <- Re(B[length(B)]); + B <- B[-length(B)]; + } res <- estimator(B, print_level); B <- res$B; CFValue <- res$CFValue; @@ -1108,7 +1116,6 @@ BICc.clm <- function(object, ...){ #' covariance matrix for the complex error. If \code{NULL} then will return value based #' on the loss used in the estimation: OLS -> "conjugate", CLS -> "direct", likelihood -> #' "matrix". -#' @param ... Other parameters passed to internal functions. #' @importFrom stats sigma #' @export sigma.clm <- function(object, type=NULL, ...){ diff --git a/man/clm.Rd b/man/clm.Rd index 7354580..41ea998 100644 --- a/man/clm.Rd +++ b/man/clm.Rd @@ -78,7 +78,39 @@ is estimated.} the data has variability and whether the regressors are correlated. Might cause trouble, especially in cases of multicollinearity.} -\item{...}{Other parameters passed to internal functions.} +\item{...}{additional parameters to pass to distribution functions. This +includes: +\itemize{ +\item \code{FI=TRUE} will make the function also produce Fisher Information +matrix, which then can be used to calculated variances of smoothing parameters +and initial states of the model. This is used in the \link[stats]{vcov} method; +\item \code{shape} - the value of the shape parameter of \code{distribution="dcgnorm"}. If +provided, the shape is not estimated; +} + +You can also pass parameters to the optimiser: +\enumerate{ +\item \code{B} - the vector of starting values of parameters for the optimiser, +should correspond to the explanatory variables. If formula for scale was provided, +the parameters for that part should follow the parameters for location; +for \code{distribution="dcgnorm"} with the estimated shape, the starting value of the shape +can be added as the last element of \code{B}; +\item \code{algorithm} - the algorithm to use in optimisation +(\code{"NLOPT_LN_SBPLX"} by default); +\item \code{maxeval} - maximum number of evaluations to carry out. Default is 40 per +estimated parameter. In case of LASSO / RIDGE the default is 80 per estimated parameter; +\item \code{maxtime} - stop, when the optimisation time (in seconds) exceeds this; +\item \code{xtol_rel} - the precision of the optimiser (the default is 1E-6); +\item \code{xtol_abs} - the absolute precision of the optimiser (the default is 1E-8); +\item \code{ftol_rel} - the stopping criterion in case of the relative change in the loss +function (the default is 1E-4); +\item \code{ftol_abs} - the stopping criterion in case of the absolute change in the loss +function (the default is 0 - not used); +\item \code{print_level} - the level of output for the optimiser (0 by default). +If equal to 41, then the detailed results of the optimisation are returned. +} +You can read more about these parameters by running the function +\link[nloptr]{nloptr.print.options}.} \item{object}{Object of class "clm" estimated via \code{clm()} function.} From 9d1eff32bcc135b037ffe8ca465cbaa06d8e548f Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 17:19:51 +0000 Subject: [PATCH 09/16] Fix cARIMA without intercept: MA sign, forecasts, AR checks, empty models - fitter() negated the MA parameters together with the ARI polynomial when there were no regressors (likelihood correct, reported coefficients wrong). - predict.clm() used parameters[1:0] (first element) and a column of ones when there was no intercept. - iModelDesign() used x[,-FALSE], selecting no columns. - estimator() now handles models with nothing to estimate. Adds regression tests. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NEWS | 9 ++++- R/clm.R | 33 +++++++++++++----- tests/testthat/test-arima-nointercept.R | 46 +++++++++++++++++++++++++ 3 files changed, 78 insertions(+), 10 deletions(-) create mode 100644 tests/testthat/test-arima-nointercept.R diff --git a/NEWS b/NEWS index 0f70021..83801eb 100644 --- a/NEWS +++ b/NEWS @@ -3,11 +3,18 @@ complex v1.0.3 (Release data: 2026-09-25) Changes: * dcgnorm() and rcgnorm() - density and random generation for the complex generalised normal distribution (complex generalised Gaussian), with scale, pseudo-scale and shape. -* clm() has a new argument distribution. distribution="dcgnorm" uses the complex generalised normal likelihood, with the shape estimated (or provided via shape=...). The scale and the pseudo-scale are concentrated out: the circularity coefficient comes from the moments of the residuals and the scale is the ML estimate given it (approximate ML when shape is not 1). Standard errors are for the parameters given the shape. +* clm() has a new argument distribution. distribution="dcgnorm" uses the complex generalised normal likelihood, with the shape estimated (or provided via shape=..., with the starting value as the last element of B). The scale and the pseudo-scale are concentrated out: the circularity coefficient comes from the moments of the residuals and the scale is the ML estimate given it (approximate ML when shape is not 1). Standard errors are for the parameters given the shape. * AIC() and BIC() of clm() now count all real parameters of the model (both parts of the complex variable), as done in multivariate models. nparam() and the degrees of freedom are still per series. Before, the penalty was half of what it should be. Information criteria values change, and so can the orders selected by them. * AICc() and BICc() methods for clm() use the multivariate small-sample correction (Bedrick & Tsai, 1994) with two series, as in the legion package. The correction assumes unrestricted multivariate regression, so it is approximate for complex models. * Added tests (testthat). +Bugfixes: +* clm() without intercept (e.g. y~-1) and with MA terms flipped the sign of the MA parameters (the likelihood was correct, the reported parameters and forecasts were not). +* predict.clm() for cARIMA without intercept used the first parameter as an intercept (1:0 indexing). +* clm() without intercept and with AR terms failed in the check of variability of differenced data. +* clm() with nothing to estimate (e.g. cARIMA(0,1,0) without intercept) failed in the optimiser. +* The documentation of the ellipsis in clm() was hidden by the one of sigma.clm(). + complex v1.0.2 (Release data: 2026-02-03) ============== diff --git a/R/clm.R b/R/clm.R index 1cec842..8ef4f1d 100644 --- a/R/clm.R +++ b/R/clm.R @@ -216,7 +216,7 @@ clm <- function(formula, data, subset, na.action, # If the vector B is not complex then it is estimated in nloptr. Make it complex if(!is.complex(B)){ nVariables <- length(B); - B <- complex(real=B[1:(nVariables/2)],imaginary=B[(nVariables/2+1):nVariables]); + B <- complex(real=B[seq_len(nVariables/2)],imaginary=B[nVariables/2+seq_len(nVariables/2)]); } # If there is ARIMA, then calculate polynomials @@ -237,7 +237,8 @@ clm <- function(formula, data, subset, na.action, B <- c(B[1:nVariablesExo], -polyprodcomplex(poly2,poly1)[-1], BMA); } else{ - B <- -c(polyprodcomplex(poly2,poly1)[-1], BMA); + # Only the ARI polynomial changes sign, not the MA parameters + B <- c(-polyprodcomplex(poly2,poly1)[-1], BMA); } } @@ -355,6 +356,11 @@ clm <- function(formula, data, subset, na.action, maxeval <- maxeval + 40; } + # Nothing to estimate (e.g. cARIMA(0,1,0) without intercept): just evaluate the loss + if(length(BReal)==0){ + return(list(B=B, CFValue=CF(B, loss, y, matrixXreg), shape=shape)); + } + # Although this is not needed in case of distribution="dnorm", we do that in a way, for the code consistency purposes res <- nloptr(BReal, CF, opts=list(algorithm=algorithm, xtol_rel=xtol_rel, maxeval=maxeval, print_level=print_level, @@ -366,7 +372,7 @@ clm <- function(formula, data, subset, na.action, if(shapeEstimate){ shapeValue <- exp(BReal[nVariables+1]); } - B[] <- complex(real=BReal[1:(nVariables/2)],imaginary=BReal[(nVariables/2+1):nVariables]); + B[] <- complex(real=BReal[seq_len(nVariables/2)],imaginary=BReal[nVariables/2+seq_len(nVariables/2)]); nVariables <- length(B); CFValue <- res$objective; @@ -822,9 +828,14 @@ clm <- function(formula, data, subset, na.action, matrixXregForDiffs <- matrixXregForDiffs[-c(1:iOrder),,drop=FALSE]; # Check variability in the new data. Have we removed important observations? - noVariability <- apply(matrixXregForDiffs[,-interceptIsNeeded,drop=FALSE]== - matrix(matrixXregForDiffs[1,-interceptIsNeeded], - nrow(matrixXregForDiffs),ncol(matrixXregForDiffs)-interceptIsNeeded, + # Columns to check: all but the intercept (x[,-FALSE] would select no columns) + columnsToCheck <- seq_len(ncol(matrixXregForDiffs)); + if(interceptIsNeeded){ + columnsToCheck <- columnsToCheck[-1]; + } + noVariability <- apply(matrixXregForDiffs[,columnsToCheck,drop=FALSE]== + matrix(matrixXregForDiffs[1,columnsToCheck], + nrow(matrixXregForDiffs),length(columnsToCheck), byrow=TRUE), 2,all); if(any(noVariability)){ @@ -832,7 +843,7 @@ clm <- function(formula, data, subset, na.action, "This might mean that all the variability for them happened ", "in the very beginning of the series. We'll try to fix this, but the model might fail.", call.=FALSE); - matrixXregForDiffs[1,which(noVariability)+1] <- rnorm(sum(noVariability)); + matrixXregForDiffs[1,columnsToCheck[noVariability]] <- rnorm(sum(noVariability)); } return(matrixXregForDiffs) @@ -1528,7 +1539,8 @@ predict.clm <- function(object, newdata=NULL, interval=c("none", "confidence", " # Split the parameters into normal and polynomial (for ARI) if(arOrderUsed || maOrderUsed){ - parameters <- parameters[1:nParametersExo]; + # seq_len() rather than 1:n, which would select the first element when n=0 + parameters <- parameters[seq_len(nParametersExo)]; } parametersNames <- names(parameters); @@ -1622,10 +1634,13 @@ predict.clm <- function(object, newdata=NULL, interval=c("none", "confidence", " matrixOfxreg <- matrixOfxreg[,parametersNames,drop=FALSE]; } else{ - matrixOfxreg <- matrix(1, h, 1); if(interceptIsNeeded){ + matrixOfxreg <- matrix(1, h, 1); colnames(matrixOfxreg) <- "(Intercept)"; } + else{ + matrixOfxreg <- matrix(0, h, 0); + } } } diff --git a/tests/testthat/test-arima-nointercept.R b/tests/testthat/test-arima-nointercept.R new file mode 100644 index 0000000..2226a54 --- /dev/null +++ b/tests/testthat/test-arima-nointercept.R @@ -0,0 +1,46 @@ +# cARIMA without intercept (formula y~-1): coefficients, residuals and forecasts +set.seed(3) +obs <- 3000 +e <- rcnorm(obs, 0, sigma2=1, varsigma2=0.3) +dy <- e +for(t in 2:obs){ + dy[t] <- (0.5+0.2i)*dy[t-1] + e[t] + (0.3-0.1i)*e[t-1] +} +y <- cumsum(dy) +data <- data.frame(y=y) + +test_that("cARIMA(1,1,1) without intercept estimates the parameters with correct signs", { + model <- suppressWarnings(clm(y~-1, data, orders=c(1,1,1))) + expect_equal(unname(coef(model)), c(0.5+0.2i, 0.3-0.1i), tolerance=0.05) + # Residuals follow the model recursion + b <- coef(model) + residualsManual <- complex(obs) + for(t in 3:obs){ + residualsManual[t] <- (y[t]-y[t-1]) - b[1]*(y[t-1]-y[t-2]) - b[2]*residualsManual[t-1] + } + expect_equal(unname(residuals(model)[21:obs]), residualsManual[21:obs], tolerance=1e-6) + # Forecasts follow the model recursion + yExtended <- y + residualsExtended <- residuals(model) + for(i in 1:5){ + n <- length(yExtended) + yExtended <- c(yExtended, yExtended[n] + b[1]*(yExtended[n]-yExtended[n-1]) + + b[2]*residualsExtended[n]) + residualsExtended <- c(residualsExtended, 0) + } + expect_equal(as.vector(predict(model, newdata=matrix(NA, 5, 1))$mean), unname(tail(yExtended, 5))) +}) + +test_that("cARIMA(0,1,2) without intercept is not affected by the sign of MA", { + model <- suppressWarnings(clm(y~-1, data, orders=c(0,1,2))) + modelIntercept <- clm(y~1, data, orders=c(0,1,2)) + expect_equal(unname(coef(model)), unname(coef(modelIntercept)[-1]), tolerance=0.02) +}) + +test_that("cARIMA(0,1,0) without intercept has nothing to estimate and forecasts the last value", { + for(distribution in c("dcnorm", "dcgnorm")){ + model <- suppressWarnings(clm(y~-1, data, orders=c(0,1,0), distribution=distribution)) + expect_true(is.finite(logLik(model))) + expect_equal(as.vector(predict(model, newdata=matrix(NA, 3, 1))$mean), rep(y[obs], 3)) + } +}) From c0489488d01365d607829ce88316d431cce1ac9d Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Fri, 25 Sep 2026 19:52:43 +0000 Subject: [PATCH 10/16] Add carima(): pure complex ARIMA wrapper of clm() with msarima-like interface Arguments y, orders (vector or list), constant (FALSE/TRUE/value), arma, model, loss, distribution, h, holdout, silent and ... passed to clm(). forecast.carima(), print.carima(), vcov.carima() (suppresses the no-intercept warning of the refit). Tests added. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NAMESPACE | 6 + NEWS | 1 + R/carima.R | 247 +++++++++++++++++++++++++++++++++++ man/carima.Rd | 103 +++++++++++++++ tests/testthat/test-carima.R | 49 +++++++ 5 files changed, 406 insertions(+) create mode 100644 R/carima.R create mode 100644 man/carima.Rd create mode 100644 tests/testthat/test-carima.R diff --git a/NAMESPACE b/NAMESPACE index 896a87a..744d1cf 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ S3method(AICc,clm) S3method(BICc,clm) S3method(actuals,clm) S3method(confint,clm) +S3method(forecast,carima) S3method(logLik,clm) S3method(nobs,clm) S3method(nparam,clm) @@ -12,11 +13,14 @@ S3method(plot,clm) S3method(plot,predict.clm) S3method(predict,clm) S3method(print,cacf) +S3method(print,carima) S3method(print,summary.clm) S3method(sigma,clm) S3method(summary,clm) +S3method(vcov,carima) S3method(vcov,clm) export(cacf) +export(carima) export(ccor) export(ccov) export(ccov2cor) @@ -61,6 +65,7 @@ importFrom(greybox, AICc, BICc, actuals, + forecast, measures, nparam, xregExpander @@ -93,6 +98,7 @@ importFrom(stats, model.matrix, nobs, pacf, + predict, qt, resid, residuals, diff --git a/NEWS b/NEWS index 83801eb..340ef78 100644 --- a/NEWS +++ b/NEWS @@ -2,6 +2,7 @@ complex v1.0.3 (Release data: 2026-09-25) ============== Changes: +* carima() - pure complex ARIMA (no regressors), a wrapper of clm() with an interface similar to smooth::msarima(): orders as c(p,d,q) or list(ar,i,ma), constant (FALSE by default, TRUE or a value), arma (fixed parameters), model (reuse of an estimated model), h, holdout (with accuracy), silent. Non-seasonal only. forecast() method for carima. * dcgnorm() and rcgnorm() - density and random generation for the complex generalised normal distribution (complex generalised Gaussian), with scale, pseudo-scale and shape. * clm() has a new argument distribution. distribution="dcgnorm" uses the complex generalised normal likelihood, with the shape estimated (or provided via shape=..., with the starting value as the last element of B). The scale and the pseudo-scale are concentrated out: the circularity coefficient comes from the moments of the residuals and the scale is the ML estimate given it (approximate ML when shape is not 1). Standard errors are for the parameters given the shape. * AIC() and BIC() of clm() now count all real parameters of the model (both parts of the complex variable), as done in multivariate models. nparam() and the degrees of freedom are still per series. Before, the penalty was half of what it should be. Information criteria values change, and so can the orders selected by them. diff --git a/R/carima.R b/R/carima.R new file mode 100644 index 0000000..dd3b28a --- /dev/null +++ b/R/carima.R @@ -0,0 +1,247 @@ +#' Complex ARIMA +#' +#' Function constructs complex ARIMA(p,d,q) model for a complex-valued time series, +#' without explanatory variables. It is a wrapper of \link[complex]{clm}. +#' +#' The model is +#' \deqn{(1-\phi_1 B-\dots-\phi_p B^p)(1-B)^d y_t = c + (1+\theta_1 B+\dots+\theta_q B^q)\epsilon_t ,} +#' where all the parameters are complex and \eqn{c} is the constant (included only when +#' \code{constant} is not \code{FALSE}). With \code{d>0} the constant acts as a drift. +#' +#' Only the non-seasonal model is supported. Parameters are estimated via \link[complex]{clm}, +#' so its restrictions apply: the first lags are extrapolated and the first MA errors are set to +#' zero (there is no \code{initial} and no \code{bounds}). Fixed parameters (\code{arma}, numeric +#' \code{constant}) can only be provided for all of the parameters of the model at once. +#' +#' @template author +#' @template keywords +#' +#' @param y vector or ts object with the complex-valued time series. +#' @param orders orders of the model: either a vector \code{c(p,d,q)} or a list +#' \code{list(ar=p, i=d, ma=q)} (elements that are not provided are set to zero). +#' @param constant if \code{TRUE}, the constant term is included in the model and estimated. +#' Can also be a complex number (the value of the constant), in which case all the other +#' parameters need to be provided via \code{arma}. +#' @param arma list \code{list(ar=..., ma=...)} or a vector (AR parameters followed by MA ones) +#' with the fixed complex AR/MA parameters. All of them need to be provided; the model is then +#' not estimated. +#' @param model previously estimated \code{carima} model. If provided, its orders, constant, +#' loss, distribution and parameters are used, and the model is applied to \code{y} without +#' re-estimation. +#' @param loss the loss function used in the estimation, see \link[complex]{clm}. +#' @param distribution the distribution used in the likelihood, see \link[complex]{clm}. +#' @param h the forecast horizon. If positive, the forecasts are produced and returned in +#' \code{forecast}. +#' @param holdout if \code{TRUE}, the last \code{h} observations are held out: the model is +#' estimated on the rest and its forecasts are compared with the held-out values. +#' @param silent if \code{FALSE}, the function plots the fit and the forecast at the end. +#' @param ... other parameters passed to \link[complex]{clm}, e.g. \code{shape} for +#' \code{distribution="dcgnorm"}, \code{B} or the parameters of the optimiser. +#' +#' @return An object of the class \code{c("carima","clm","greybox")}: the model returned by +#' \link[complex]{clm}, plus +#' \itemize{ +#' \item \code{orders} - the list of orders, +#' \item \code{constant} - whether the constant is included, +#' \item \code{h}, \code{holdout} - the horizon and the holdout flag, +#' \item \code{forecast} - the point forecasts (if \code{h>0}), +#' \item \code{holdoutData} - the held-out values (if \code{holdout=TRUE}), +#' \item \code{accuracy} - the error measures on the holdout: real and imaginary parts of the +#' mean error (ME_r, ME_i), MAE (mean modulus of the error) and RMSE (root mean squared modulus +#' of the error). +#' } +#' +#' @seealso \code{\link[complex]{clm}}, \code{\link[smooth]{msarima}} +#' +#' @examples +#' set.seed(41) +#' e <- rcnorm(200, 0, sigma2=1, varsigma2=0.3+0.2i) +#' y <- cumsum(e + c(0, 0.5*e[-200])) +#' model <- carima(y, orders=c(0,1,1), h=10, holdout=TRUE) +#' summary(model) +#' model$accuracy +#' forecast(model, h=5) +#' +#' @rdname carima +#' @export +carima <- function(y, orders=list(ar=0, i=0, ma=0), constant=FALSE, arma=NULL, model=NULL, + loss=c("likelihood","OLS","CLS","MSE","MAE","HAM"), + distribution=c("dcnorm","dcgnorm"), + h=0, holdout=FALSE, silent=TRUE, ...){ + cl <- match.call(); + ellipsis <- list(...); + + # Reuse the previously estimated model + if(!is.null(model)){ + if(!inherits(model, "carima")){ + stop("model should be estimated via carima().", call.=FALSE); + } + orders <- model$orders; + constant <- model$constant; + loss <- model$loss; + distribution <- model$distribution; + parameters <- coef(model); + if(distribution=="dcgnorm"){ + ellipsis$shape <- model$other$shape; + } + } + else{ + loss <- match.arg(loss); + distribution <- match.arg(distribution); + parameters <- NULL; + } + + orders <- carimaOrders(orders); + ordersVector <- c(orders$ar, orders$i, orders$ma); + + # Constant: FALSE, TRUE (estimate) or a value + if(is.logical(constant)){ + constantValue <- NULL; + } + else{ + constantValue <- as.complex(constant); + constant <- TRUE; + } + + # Fixed parameters: all or nothing + if(is.null(model) && (!is.null(arma) || !is.null(constantValue))){ + if(is.list(arma)){ + arma <- c(arma$ar, arma$ma); + } + arma <- as.complex(arma); + if(length(arma)!=orders$ar+orders$ma){ + stop("arma should contain all ", orders$ar+orders$ma, " AR/MA parameters.", call.=FALSE); + } + if(constant && is.null(constantValue)){ + stop("The constant needs to be provided as a value when arma is provided, ", + "because clm() cannot estimate a part of the parameters.", call.=FALSE); + } + parameters <- c(constantValue, arma); + } + + # Holdout + y <- as.vector(y); + if(!is.complex(y)){ + y <- as.complex(y); + } + obsAll <- length(y); + if(holdout){ + if(h<=0){ + stop("holdout=TRUE requires h>0.", call.=FALSE); + } + holdoutData <- y[obsAll-h+seq_len(h)]; + y <- y[seq_len(obsAll-h)]; + } + + # Estimate via clm + formula <- if(constant) y~1 else y~-1; + clmArgs <- c(list(formula=formula, data=data.frame(y=y), loss=loss, distribution=distribution, + orders=ordersVector, parameters=parameters), + ellipsis); + # clm warns about the missing intercept, which is intended in carima() + ourModel <- suppressWarnings(do.call(clm, clmArgs)); + if(!is.null(parameters)){ + names(ourModel$coefficients) <- names(ourModel$B) <- carimaNames(orders, constant); + } + + # The clm call is kept in ourModel$call, because vcov.clm() re-estimates the model from it + ourModel$carimaCall <- cl; + ourModel$orders <- orders; + ourModel$constant <- constant; + ourModel$h <- h; + ourModel$holdout <- holdout; + ourModel$model <- paste0("cARIMA(", paste0(ordersVector, collapse=","), ")", + ifelse(constant, " with constant", "")); + class(ourModel) <- c("carima", class(ourModel)); + + if(h>0){ + ourModel$forecast <- as.vector(predict(ourModel, newdata=matrix(NA, h, 1))$mean); + if(holdout){ + ourModel$holdoutData <- holdoutData; + errors <- holdoutData - ourModel$forecast; + ourModel$accuracy <- c(ME_r=Re(mean(errors)), ME_i=Im(mean(errors)), + MAE=mean(Mod(errors)), RMSE=sqrt(mean(Mod(errors)^2))); + } + } + + if(!silent){ + if(h>0){ + plot(forecast(ourModel, h=h)); + } + else{ + plot(ourModel, 7); + } + } + + return(ourModel); +} + +# Orders as a list(ar, i, ma) with single non-negative values +carimaOrders <- function(orders){ + if(is.list(orders)){ + orders <- list(ar=if(is.null(orders$ar)) 0 else orders$ar, + i=if(is.null(orders$i)) 0 else orders$i, + ma=if(is.null(orders$ma)) 0 else orders$ma); + } + else{ + if(length(orders)!=3){ + stop("orders should be a vector c(p,d,q) or a list(ar=p, i=d, ma=q).", call.=FALSE); + } + orders <- list(ar=orders[1], i=orders[2], ma=orders[3]); + } + if(any(sapply(orders, length)!=1)){ + stop("Only non-seasonal cARIMA is supported: each order should be a single number.", + call.=FALSE); + } + if(any(unlist(orders)<0)){ + stop("Orders should be non-negative.", call.=FALSE); + } + return(orders); +} + +# Parameter names as clm produces them +carimaNames <- function(orders, constant){ + c(if(constant) "(Intercept)", + if(orders$ar>0) paste0("yLag", seq_len(orders$ar)), + if(orders$ma>0) paste0("eLag", seq_len(orders$ma))); +} + +#' @param object the model estimated via \code{carima()}. +#' @param level the confidence level for the prediction interval. +#' @param interval the type of interval, see \link[complex]{clm}. +#' @rdname carima +#' @importFrom greybox forecast +#' @importFrom stats predict +#' @export +forecast.carima <- function(object, h=10, interval=c("none", "confidence", "prediction"), + level=0.95, ...){ + interval <- match.arg(interval); + return(predict(object, newdata=matrix(NA, h, 1), interval=interval, level=level, ...)); +} + +#' @export +print.carima <- function(x, ...){ + cat(x$model, "estimated via carima()\n"); + cat("Loss function used in estimation:", x$loss, "\n"); + if(x$distribution=="dcgnorm"){ + cat("Distribution: complex generalised normal, shape =", round(x$other$shape, 4), "\n"); + } + cat("Coefficients:\n"); + print(coef(x)); + if(!is.na(logLik(x))){ + ICs <- c(AIC=AIC(x), AICc=AICc(x), BIC=BIC(x), BICc=BICc(x)); + cat("Information criteria:\n"); + print(round(ICs, 4)); + } + if(!is.null(x$accuracy)){ + cat("Holdout accuracy (h=", x$h, "):\n", sep=""); + print(x$accuracy); + } +} + +# vcov.clm() re-estimates the model via clm(), which warns about the missing intercept. +# In carima() this is intended, so the warning is suppressed. +#' @export +vcov.carima <- function(object, ...){ + return(suppressWarnings(NextMethod())); +} diff --git a/man/carima.Rd b/man/carima.Rd new file mode 100644 index 0000000..ffbf092 --- /dev/null +++ b/man/carima.Rd @@ -0,0 +1,103 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/carima.R +\name{carima} +\alias{carima} +\alias{forecast.carima} +\title{Complex ARIMA} +\usage{ +carima(y, orders = list(ar = 0, i = 0, ma = 0), constant = FALSE, + arma = NULL, model = NULL, loss = c("likelihood", "OLS", "CLS", "MSE", + "MAE", "HAM"), distribution = c("dcnorm", "dcgnorm"), h = 0, + holdout = FALSE, silent = TRUE, ...) + +\method{forecast}{carima}(object, h = 10, interval = c("none", + "confidence", "prediction"), level = 0.95, ...) +} +\arguments{ +\item{y}{vector or ts object with the complex-valued time series.} + +\item{orders}{orders of the model: either a vector \code{c(p,d,q)} or a list +\code{list(ar=p, i=d, ma=q)} (elements that are not provided are set to zero).} + +\item{constant}{if \code{TRUE}, the constant term is included in the model and estimated. +Can also be a complex number (the value of the constant), in which case all the other +parameters need to be provided via \code{arma}.} + +\item{arma}{list \code{list(ar=..., ma=...)} or a vector (AR parameters followed by MA ones) +with the fixed complex AR/MA parameters. All of them need to be provided; the model is then +not estimated.} + +\item{model}{previously estimated \code{carima} model. If provided, its orders, constant, +loss, distribution and parameters are used, and the model is applied to \code{y} without +re-estimation.} + +\item{loss}{the loss function used in the estimation, see \link[complex]{clm}.} + +\item{distribution}{the distribution used in the likelihood, see \link[complex]{clm}.} + +\item{h}{the forecast horizon. If positive, the forecasts are produced and returned in +\code{forecast}.} + +\item{holdout}{if \code{TRUE}, the last \code{h} observations are held out: the model is +estimated on the rest and its forecasts are compared with the held-out values.} + +\item{silent}{if \code{FALSE}, the function plots the fit and the forecast at the end.} + +\item{...}{other parameters passed to \link[complex]{clm}, e.g. \code{shape} for +\code{distribution="dcgnorm"}, \code{B} or the parameters of the optimiser.} + +\item{object}{the model estimated via \code{carima()}.} + +\item{interval}{the type of interval, see \link[complex]{clm}.} + +\item{level}{the confidence level for the prediction interval.} +} +\value{ +An object of the class \code{c("carima","clm","greybox")}: the model returned by +\link[complex]{clm}, plus +\itemize{ +\item \code{orders} - the list of orders, +\item \code{constant} - whether the constant is included, +\item \code{h}, \code{holdout} - the horizon and the holdout flag, +\item \code{forecast} - the point forecasts (if \code{h>0}), +\item \code{holdoutData} - the held-out values (if \code{holdout=TRUE}), +\item \code{accuracy} - the error measures on the holdout: real and imaginary parts of the +mean error (ME_r, ME_i), MAE (mean modulus of the error) and RMSE (root mean squared modulus +of the error). +} +} +\description{ +Function constructs complex ARIMA(p,d,q) model for a complex-valued time series, +without explanatory variables. It is a wrapper of \link[complex]{clm}. +} +\details{ +The model is +\deqn{(1-\phi_1 B-\dots-\phi_p B^p)(1-B)^d y_t = c + (1+\theta_1 B+\dots+\theta_q B^q)\epsilon_t ,} +where all the parameters are complex and \eqn{c} is the constant (included only when +\code{constant} is not \code{FALSE}). With \code{d>0} the constant acts as a drift. + +Only the non-seasonal model is supported. Parameters are estimated via \link[complex]{clm}, +so its restrictions apply: the first lags are extrapolated and the first MA errors are set to +zero (there is no \code{initial} and no \code{bounds}). Fixed parameters (\code{arma}, numeric +\code{constant}) can only be provided for all of the parameters of the model at once. +} +\examples{ +set.seed(41) +e <- rcnorm(200, 0, sigma2=1, varsigma2=0.3+0.2i) +y <- cumsum(e + c(0, 0.5*e[-200])) +model <- carima(y, orders=c(0,1,1), h=10, holdout=TRUE) +summary(model) +model$accuracy +forecast(model, h=5) + +} +\seealso{ +\code{\link[complex]{clm}}, \code{\link[smooth]{msarima}} +} +\author{ +Ivan Svetunkov, \email{ivan@svetunkov.com} +} +\keyword{models} +\keyword{nonlinear} +\keyword{regression} +\keyword{ts} diff --git a/tests/testthat/test-carima.R b/tests/testthat/test-carima.R new file mode 100644 index 0000000..db11f56 --- /dev/null +++ b/tests/testthat/test-carima.R @@ -0,0 +1,49 @@ +# carima(): wrapper of clm() for pure complex ARIMA +set.seed(41) +obs <- 300 +e <- rcnorm(obs, 0, sigma2=1, varsigma2=0.3+0.2i) +y <- cumsum(e + c(0, (0.5-0.2i)*e[-obs])) + +test_that("carima() gives the same model as clm() without intercept", { + model <- carima(y, orders=c(0,1,1)) + modelClm <- suppressWarnings(clm(y~-1, data.frame(y=y), orders=c(0,1,1))) + expect_s3_class(model, "carima") + expect_equal(unname(coef(model)), unname(coef(modelClm))) + expect_equal(as.numeric(logLik(model)), as.numeric(logLik(modelClm))) + expect_equal(model$orders, list(ar=0, i=1, ma=1)) +}) + +test_that("orders can be a list and the constant is included when asked", { + model <- carima(y, orders=list(i=1, ma=1), constant=TRUE) + modelClm <- clm(y~1, data.frame(y=y), orders=c(0,1,1)) + expect_equal(unname(coef(model)), unname(coef(modelClm))) + expect_error(carima(y, orders=list(ar=c(1,1))), "non-seasonal") +}) + +test_that("holdout and forecast work", { + model <- carima(y, orders=c(0,1,1), h=10, holdout=TRUE) + expect_equal(nobs(model), obs - 10) + expect_length(model$forecast, 10) + expect_equal(model$holdoutData, y[obs-9:0]) + expect_equal(unname(model$accuracy["RMSE"]), + sqrt(mean(Mod(model$holdoutData - model$forecast)^2))) + expect_equal(as.vector(forecast(model, h=10)$mean), model$forecast) + expect_no_warning(summary(model)) +}) + +test_that("model and arma reuse parameters without estimation", { + model <- carima(y, orders=c(0,1,1)) + modelReused <- carima(y[1:250], model=model) + expect_equal(coef(modelReused), coef(model)) + modelFixed <- carima(y, orders=c(0,1,1), arma=list(ma=0.5-0.2i)) + expect_equal(unname(coef(modelFixed)), 0.5-0.2i) + expect_error(carima(y, orders=c(0,1,1), constant=TRUE, arma=list(ma=0.5)), "constant") + expect_error(carima(y, orders=c(1,1,1), arma=list(ma=0.5)), "all") +}) + +test_that("carima() works with the complex generalised normal", { + model <- carima(y, orders=c(0,1,1), distribution="dcgnorm") + expect_true(model$other$shapeEstimated) + modelReused <- carima(y[1:250], model=model) + expect_equal(modelReused$other$shape, model$other$shape) +}) From 2f037f7a3f88107da970c93ce720d5265dfecbeb Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Sat, 26 Sep 2026 12:29:57 +0000 Subject: [PATCH 11/16] Add auto.carima() and stationarity/invertibility penalty in clm() auto.carima(): AR/MA orders selected per d via neighbourhood search (auto.ssarima principle of selecting ARMA within each d), constant tested at the end for the best model of each d with the search continued from there; search="full" for all combinations. ICs on the common sample without the first ar+i+ma observations (models estimated on the full sample). Warm starts from the closest estimated model. clm(): loss returns 1E+100/min|root| if the complex AR or MA polynomial has roots inside the unit circle. Information criteria formulas moved to clmIC(). carima() forecasting/holdout moved to carimaFinalise(). Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NAMESPACE | 1 + NEWS | 2 + R/autoCarima.R | 298 ++++++++++++++++++++++++++++++ R/carima.R | 11 +- R/clm.R | 74 +++++--- man/auto.carima.Rd | 104 +++++++++++ tests/testthat/test-auto-carima.R | 57 ++++++ 7 files changed, 520 insertions(+), 27 deletions(-) create mode 100644 R/autoCarima.R create mode 100644 man/auto.carima.Rd create mode 100644 tests/testthat/test-auto-carima.R diff --git a/NAMESPACE b/NAMESPACE index 744d1cf..88cb2a7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,6 +19,7 @@ S3method(sigma,clm) S3method(summary,clm) S3method(vcov,carima) S3method(vcov,clm) +export(auto.carima) export(cacf) export(carima) export(ccor) diff --git a/NEWS b/NEWS index 340ef78..8b99125 100644 --- a/NEWS +++ b/NEWS @@ -3,6 +3,8 @@ complex v1.0.3 (Release data: 2026-09-25) Changes: * carima() - pure complex ARIMA (no regressors), a wrapper of clm() with an interface similar to smooth::msarima(): orders as c(p,d,q) or list(ar,i,ma), constant (FALSE by default, TRUE or a value), arma (fixed parameters), model (reuse of an estimated model), h, holdout (with accuracy), silent. Non-seasonal only. forecast() method for carima. +* auto.carima() - automatic selection of the orders of complex ARIMA. Following auto.ssarima(), AR and MA orders are selected for each order of differencing via an IC-driven neighbourhood search (Hyndman & Khandakar, 2008), and the constant is tested at the end for the best model of each d. search="full" fits all the combinations. The information criteria are calculated on the common sample without the first ar+i+ma observations, so that models with different d are comparable. Each model is estimated from the default starting values and from the closest estimated model. +* clm() now penalises non-stationary AR and non-invertible MA parameters in the loss (all roots of the complex AR/MA polynomials should lie outside the unit circle); the loss becomes 1E+100 divided by the smallest modulus of the roots. * dcgnorm() and rcgnorm() - density and random generation for the complex generalised normal distribution (complex generalised Gaussian), with scale, pseudo-scale and shape. * clm() has a new argument distribution. distribution="dcgnorm" uses the complex generalised normal likelihood, with the shape estimated (or provided via shape=..., with the starting value as the last element of B). The scale and the pseudo-scale are concentrated out: the circularity coefficient comes from the moments of the residuals and the scale is the ML estimate given it (approximate ML when shape is not 1). Standard errors are for the parameters given the shape. * AIC() and BIC() of clm() now count all real parameters of the model (both parts of the complex variable), as done in multivariate models. nparam() and the degrees of freedom are still per series. Before, the penalty was half of what it should be. Information criteria values change, and so can the orders selected by them. diff --git a/R/autoCarima.R b/R/autoCarima.R new file mode 100644 index 0000000..3556ce7 --- /dev/null +++ b/R/autoCarima.R @@ -0,0 +1,298 @@ +#' Automatic Complex ARIMA +#' +#' Function selects the orders of complex ARIMA (see \link[complex]{carima}) based on an +#' information criterion. +#' +#' The selection follows the principles of \link[smooth]{auto.ssarima}: the AR and MA orders are +#' selected for each order of differencing, and the best model across all of them is chosen. +#' With \code{search="stepwise"}: +#' \enumerate{ +#' \item For each \eqn{d} from 0 to \code{orders$i}, the AR and MA orders are selected via a +#' neighbourhood search (similar to Hyndman & Khandakar, 2008): starting from cARIMA(0,d,0), +#' all models with p and q different by at most one are fitted, and the search moves to the best +#' of them, as long as the information criterion improves. If \code{constant=NULL}, the constant +#' is included for \eqn{d=0} and excluded otherwise; +#' \item If \code{constant=NULL}, the best model of each \eqn{d} is re-estimated with the opposite +#' option for the constant and the neighbourhood search continues from it. The best of all the +#' estimated models is returned. +#' } +#' With \code{search="full"}, all combinations of p, d, q and the constant are fitted +#' (except for cARIMA(p,0,q) without constant when \code{constant=NULL}). +#' +#' Models with different orders of differencing and AR use different numbers of initial values, +#' which \link[complex]{clm} constructs from the data. To make the information criteria +#' comparable, all models are estimated on the whole sample, but their criteria are calculated on +#' the common sample without the first \eqn{k} observations, where \eqn{k} is the sum of the +#' maximum orders. The log-likelihood on that sample is calculated at the estimated parameters, +#' with the scale re-estimated. +#' +#' Every model is estimated from the default starting values and from the parameters of the +#' closest model already estimated (with zeros for the new parameters), and the one with the +#' higher likelihood is kept. +#' +#' @template author +#' @template keywords +#' +#' @param y vector or ts object with the complex-valued time series. +#' @param orders list of maximum orders \code{list(ar=p, i=d, ma=q)} (or vector \code{c(p,d,q)}). +#' @param constant if \code{NULL}, the function checks whether the constant is needed. +#' Otherwise \code{TRUE} or \code{FALSE}. +#' @param ic the information criterion used in the selection. +#' @param search \code{"stepwise"} for the neighbourhood search or \code{"full"} for all the +#' combinations of orders. +#' @param distribution the distribution used in the likelihood, see \link[complex]{clm}. +#' @param h the forecast horizon. +#' @param holdout if \code{TRUE}, the last \code{h} observations are held out. +#' @param silent if \code{FALSE}, the progress is printed and the final model is plotted. +#' @param ... other parameters passed to \link[complex]{carima} and \link[complex]{clm}. +#' +#' @return The selected model of the class \code{carima} (see \link[complex]{carima}), with +#' additional elements: +#' \itemize{ +#' \item \code{ICs} - data frame with the orders, the constant and the information criterion +#' (on the common sample) of all the models estimated in the process, +#' \item \code{ic} - the information criterion used, +#' \item \code{icObsDropped} - the number of first observations dropped for the criteria. +#' } +#' +#' @references \itemize{ +#' \item Svetunkov, I., Boylan, J. E. (2019). State-space ARIMA for supply-chain forecasting. +#' International Journal of Production Research, 58(3), 818-827. +#' \item Hyndman, R. J., Khandakar, Y. (2008). Automatic time series forecasting: the forecast +#' package for R. Journal of Statistical Software, 27(3), 1-22. +#' } +#' +#' @seealso \code{\link[complex]{carima}}, \code{\link[smooth]{auto.ssarima}}, +#' \code{\link[smooth]{auto.msarima}} +#' +#' @examples +#' set.seed(41) +#' e <- rcnorm(300, 0, sigma2=1, varsigma2=0.3+0.2i) +#' y <- cumsum(e + c(0, (0.5-0.2i)*e[-300])) +#' model <- auto.carima(y, orders=list(ar=2, i=1, ma=2)) +#' model +#' model$ICs +#' +#' @rdname auto.carima +#' @export +auto.carima <- function(y, orders=list(ar=3, i=2, ma=3), constant=NULL, + ic=c("AICc","AIC","BIC","BICc"), search=c("stepwise","full"), + distribution=c("dcnorm","dcgnorm"), + h=0, holdout=FALSE, silent=TRUE, ...){ + cl <- match.call(); + ic <- match.arg(ic); + search <- match.arg(search); + distribution <- match.arg(distribution); + ellipsis <- list(...); + # The starting values are set by the function + ellipsis$B <- NULL; + + orders <- carimaOrders(orders); + arMax <- orders$ar; + iMax <- orders$i; + maMax <- orders$ma; + + y <- as.vector(y); + if(!is.complex(y)){ + y <- as.complex(y); + } + obsAll <- length(y); + holdoutData <- NULL; + if(holdout){ + if(h<=0){ + stop("holdout=TRUE requires h>0.", call.=FALSE); + } + holdoutData <- y[obsAll-h+seq_len(h)]; + y <- y[seq_len(obsAll-h)]; + } + + # Observations dropped for the calculation of the information criteria + obsDropped <- arMax + iMax + maMax; + if(length(y)-obsDropped <= 2*(arMax+maMax+1)+3){ + stop("The sample is too small for the maximum orders provided.", call.=FALSE); + } + + # The information criterion on the common sample + icTrimmed <- function(model){ + errors <- residuals(model)[-seq_len(obsDropped)]; + obsTrimmed <- length(errors); + if(distribution=="dcnorm"){ + errorsMatrix <- complex2vec(errors); + scaleMatrix <- t(errorsMatrix) %*% errorsMatrix / obsTrimmed; + llikelihood <- -obsTrimmed*(log(2*pi) + 1 + 0.5*log(det(scaleMatrix))); + } + else{ + llikelihood <- cgnormConcentrated(errors, model$other$shape)$logLik; + } + return(clmIC(llikelihood, nparam(model), obsTrimmed, + isTRUE(model$other$shapeEstimated), ic)); + } + + # Starting values from an estimated model for cARIMA(p,d,q) with the given constant + startFrom <- function(model, p, q, constant){ + coefs <- coef(model); + pFrom <- model$orders$ar; + qFrom <- model$orders$ma; + k <- as.integer(model$constant); + start <- c(if(constant) {if(model$constant) coefs[1] else 0+0i}, + coefs[k+seq_len(min(p, pFrom))], rep(0+0i, max(p-pFrom, 0)), + coefs[k+pFrom+seq_len(min(q, qFrom))], rep(0+0i, max(q-qFrom, 0))); + if(distribution=="dcgnorm" && isTRUE(model$other$shapeEstimated)){ + start <- c(start, model$other$shape); + } + return(unname(start)); + } + + # Estimated models: list of list(model, IC), named by the orders and the constant + estimated <- list(); + modelKey <- function(p, d, q, constant){ + paste(p, d, q, constant); + } + fitCandidate <- function(p, d, q, constant, startModels=list()){ + key <- modelKey(p, d, q, constant); + if(!is.null(estimated[[key]])){ + return(estimated[[key]]); + } + fit <- function(B=NULL){ + carimaArgs <- c(list(y=y, orders=c(p, d, q), constant=constant, + loss="likelihood", distribution=distribution), + if(!is.null(B)) list(B=B), ellipsis); + tryCatch(do.call(carima, carimaArgs), error=function(e) NULL); + } + fits <- c(list(fit()), + lapply(startModels, function(model) fit(startFrom(model, p, q, constant)))); + fits <- fits[!sapply(fits, is.null)]; + if(length(fits)==0){ + result <- list(model=NULL, IC=Inf); + } + else{ + model <- fits[[which.max(sapply(fits, function(x) as.numeric(logLik(x))))]]; + ICValue <- icTrimmed(model); + result <- list(model=model, IC=ifelse(is.finite(ICValue), ICValue, Inf)); + } + estimated[[key]] <<- result; + if(!silent){ + cat("cARIMA(", p, ",", d, ",", q, ")", ifelse(constant, " with constant", ""), + ": ", round(result$IC, 3), "\n", sep=""); + } + return(result); + } + + constantOptions <- if(is.null(constant)) c(FALSE, TRUE) else constant; + + if(search=="full"){ + for(d in 0:iMax){ + for(constantValue in constantOptions){ + if(d==0 && !constantValue && is.null(constant)){ + next; + } + # By level p+q, so that each model can start from the nested ones + for(level in 0:(arMax+maMax)){ + for(p in 0:arMax){ + q <- level - p; + if(q<0 || q>maMax){ + next; + } + startModels <- list(); + if(p>0){ + startModels$ar <- estimated[[modelKey(p-1, d, q, constantValue)]]$model; + } + if(q>0){ + startModels$ma <- estimated[[modelKey(p, d, q-1, constantValue)]]$model; + } + fitCandidate(p, d, q, constantValue, startModels[!sapply(startModels, is.null)]); + } + } + } + } + } + else{ + # Neighbourhood search: from the given model, move to the best of the models with p and q + # different by at most one, as long as the information criterion improves + neighbourhoodSearch <- function(best, d, constantValue){ + pBest <- best$model$orders$ar; + qBest <- best$model$orders$ma; + repeat{ + neighbours <- expand.grid(p=pBest+(-1:1), q=qBest+(-1:1)); + neighbours <- neighbours[neighbours$p>=0 & neighbours$p<=arMax & + neighbours$q>=0 & neighbours$q<=maMax & + !(neighbours$p==pBest & neighbours$q==qBest),,drop=FALSE]; + if(nrow(neighbours)==0){ + break; + } + results <- lapply(seq_len(nrow(neighbours)), function(i){ + fitCandidate(neighbours$p[i], d, neighbours$q[i], constantValue, + if(!is.null(best$model)) list(best$model) else list()); + }); + neighbourICs <- sapply(results, `[[`, "IC"); + if(min(neighbourICs) < best$IC){ + best <- results[[which.min(neighbourICs)]]; + pBest <- neighbours$p[which.min(neighbourICs)]; + qBest <- neighbours$q[which.min(neighbourICs)]; + } + else{ + break; + } + } + return(best); + } + + # Search for p and q within each d, starting from cARIMA(0,d,0). The constant is included + # for d=0 and excluded for d>0 (unless it is provided), and is tested at the end. + bestPerDiff <- vector("list", iMax+1); + for(d in 0:iMax){ + constantValue <- if(is.null(constant)) (d==0) else constant; + if(!silent){ + cat("Searching AR and MA orders with d=", d, "...\n", sep=""); + } + best <- fitCandidate(0, d, 0, constantValue); + if(!is.null(best$model)){ + best <- neighbourhoodSearch(best, d, constantValue); + } + bestPerDiff[[d+1]] <- best; + } + + # Test the constant on the best model of each d, and continue the search from there: + # e.g. a random walk with drift is otherwise missed, when the search without constant + # compensates the drift with MA terms + if(is.null(constant)){ + if(!silent){ + cat("Testing the constant...\n"); + } + for(d in 0:iMax){ + bestPerD <- bestPerDiff[[d+1]]; + if(!is.null(bestPerD$model)){ + constantValue <- !bestPerD$model$constant; + flipped <- fitCandidate(bestPerD$model$orders$ar, d, bestPerD$model$orders$ma, + constantValue, list(bestPerD$model)); + if(!is.null(flipped$model)){ + neighbourhoodSearch(flipped, d, constantValue); + } + } + } + } + } + + # The best model among all the estimated ones + ICs <- sapply(estimated, `[[`, "IC"); + if(all(!is.finite(ICs))){ + stop("None of the models could be estimated.", call.=FALSE); + } + bestModel <- estimated[[which.min(ICs)]]$model; + + ICsTable <- do.call(rbind, lapply(strsplit(names(estimated), " "), function(x){ + data.frame(p=as.integer(x[1]), d=as.integer(x[2]), q=as.integer(x[3]), + constant=as.logical(x[4])); + })); + ICsTable[[ic]] <- ICs; + ICsTable <- ICsTable[order(ICsTable[[ic]]),]; + rownames(ICsTable) <- NULL; + + bestModel$ICs <- ICsTable; + bestModel$ic <- ic; + bestModel$icObsDropped <- obsDropped; + bestModel$autoCall <- cl; + bestModel <- carimaFinalise(bestModel, h, holdout, holdoutData, silent); + + return(bestModel); +} diff --git a/R/carima.R b/R/carima.R index dd3b28a..46e4736 100644 --- a/R/carima.R +++ b/R/carima.R @@ -154,6 +154,15 @@ carima <- function(y, orders=list(ar=0, i=0, ma=0), constant=FALSE, arma=NULL, m ifelse(constant, " with constant", "")); class(ourModel) <- c("carima", class(ourModel)); + ourModel <- carimaFinalise(ourModel, h, holdout, if(holdout) holdoutData else NULL, silent); + + return(ourModel); +} + +# Forecasts, holdout accuracy and the plot for an estimated carima model +carimaFinalise <- function(ourModel, h, holdout, holdoutData, silent){ + ourModel$h <- h; + ourModel$holdout <- holdout; if(h>0){ ourModel$forecast <- as.vector(predict(ourModel, newdata=matrix(NA, h, 1))$mean); if(holdout){ @@ -163,7 +172,6 @@ carima <- function(y, orders=list(ar=0, i=0, ma=0), constant=FALSE, arma=NULL, m MAE=mean(Mod(errors)), RMSE=sqrt(mean(Mod(errors)^2))); } } - if(!silent){ if(h>0){ plot(forecast(ourModel, h=h)); @@ -172,7 +180,6 @@ carima <- function(y, orders=list(ar=0, i=0, ma=0), constant=FALSE, arma=NULL, m plot(ourModel, 7); } } - return(ourModel); } diff --git a/R/clm.R b/R/clm.R index 8ef4f1d..6f0e93a 100644 --- a/R/clm.R +++ b/R/clm.R @@ -283,6 +283,27 @@ clm <- function(formula, data, subset, na.action, shapeValue <- exp(B[length(B)]); B <- B[-length(B)]; } + + # Stationarity (AR) and invertibility (MA) conditions: all roots of the complex polynomials + # should lie outside the unit circle. The differencing polynomial is not checked. + # The check is skipped if the vector of parameters does not match the orders (this happens + # in the re-estimation done by vcov.clm, where the MA columns are passed as regressors) + if((arOrder>0 || maOrder>0) && + (ifelse(is.complex(B), length(B), length(B)/2) == nVariablesExo+arOrder+maOrder)){ + if(is.complex(B)){ + BComplex <- B; + } + else{ + BComplex <- complex(real=B[seq_len(length(B)/2)], + imaginary=B[length(B)/2+seq_len(length(B)/2)]); + } + rootsModuli <- c(if(arOrder>0) Mod(polyroot(c(1, -BComplex[nVariablesExo+seq_len(arOrder)]))), + if(maOrder>0) Mod(polyroot(c(1, BComplex[nVariablesExo+arOrder+seq_len(maOrder)])))); + if(any(rootsModuli<=1)){ + return(1/min(rootsModuli)*1E+100); + } + } + fitterReturn <- fitter(B, y, matrixXreg); if(loss=="likelihood" && distribution=="dcgnorm"){ @@ -1082,41 +1103,44 @@ logLik.clm <- function(object, ...){ } # Multivariate small-sample corrections (Bedrick & Tsai, 1994), as in legion, -# with two series (real and imaginary parts). nParamPerSeries excludes the share of -# the covariance matrix (3 real parameters, i.e. 3/2 per series). The correction is -# derived for unrestricted multivariate regression, so it is approximate for clm. +# with two series (real and imaginary parts). The correction is derived for unrestricted +# multivariate regression, so it is approximate for clm. #' @export AICc.clm <- function(object, ...){ - llikelihood <- logLik(object); - llikelihood <- llikelihood[1:length(llikelihood)]; - nSeries <- 2; - nParamAll <- 2*nparam(object); - nParamPerSeries <- nparam(object) - (nSeries+1)/2 - isTRUE(object$other$shapeEstimated)/2; - obs <- nobs(object); - if(obs - (nParamPerSeries + nSeries + 1) <= 0){ - IC <- Inf; - } - else{ - IC <- -2*llikelihood + 2*(obs*nParamAll/(obs - (nParamPerSeries + nSeries + 1))); - } - return(IC); + return(clmIC(logLik(object), nparam(object), nobs(object), + isTRUE(object$other$shapeEstimated), "AICc")); } #' @export BICc.clm <- function(object, ...){ - llikelihood <- logLik(object); - llikelihood <- llikelihood[1:length(llikelihood)]; + return(clmIC(logLik(object), nparam(object), nobs(object), + isTRUE(object$other$shapeEstimated), "BICc")); +} + +# Information criteria of clm given the log-likelihood, the number of parameters per series +# (as returned by nparam.clm) and the number of observations. +# All real parameters (2*nParam) are counted in the penalty. For AICc/BICc, the number of +# parameters per series excludes the share of the covariance matrix (3/2) and of the shape (1/2). +clmIC <- function(llikelihood, nParam, obs, shapeEstimated=FALSE, ic=c("AIC","AICc","BIC","BICc")){ + ic <- match.arg(ic); + llikelihood <- as.numeric(llikelihood); nSeries <- 2; - nParamAll <- 2*nparam(object); - nParamPerSeries <- nparam(object) - (nSeries+1)/2 - isTRUE(object$other$shapeEstimated)/2; - obs <- nobs(object); + nParamAll <- 2*nParam; + nParamPerSeries <- nParam - (nSeries+1)/2 - shapeEstimated/2; + if(ic=="AIC"){ + return(-2*llikelihood + 2*nParamAll); + } + if(ic=="BIC"){ + return(-2*llikelihood + log(obs)*nParamAll); + } if(obs - (nParamPerSeries + nSeries + 1) <= 0){ - IC <- Inf; + return(Inf); } - else{ - IC <- -2*llikelihood + log(obs)*(obs*nParamAll/(obs - (nParamPerSeries + nSeries + 1))); + penalty <- obs*nParamAll/(obs - (nParamPerSeries + nSeries + 1)); + if(ic=="AICc"){ + return(-2*llikelihood + 2*penalty); } - return(IC); + return(-2*llikelihood + log(obs)*penalty); } #' @rdname clm diff --git a/man/auto.carima.Rd b/man/auto.carima.Rd new file mode 100644 index 0000000..cf0cfae --- /dev/null +++ b/man/auto.carima.Rd @@ -0,0 +1,104 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/autoCarima.R +\name{auto.carima} +\alias{auto.carima} +\title{Automatic Complex ARIMA} +\usage{ +auto.carima(y, orders = list(ar = 3, i = 2, ma = 3), constant = NULL, + ic = c("AICc", "AIC", "BIC", "BICc"), search = c("stepwise", "full"), + distribution = c("dcnorm", "dcgnorm"), h = 0, holdout = FALSE, + silent = TRUE, ...) +} +\arguments{ +\item{y}{vector or ts object with the complex-valued time series.} + +\item{orders}{list of maximum orders \code{list(ar=p, i=d, ma=q)} (or vector \code{c(p,d,q)}).} + +\item{constant}{if \code{NULL}, the function checks whether the constant is needed. +Otherwise \code{TRUE} or \code{FALSE}.} + +\item{ic}{the information criterion used in the selection.} + +\item{search}{\code{"stepwise"} for the neighbourhood search or \code{"full"} for all the +combinations of orders.} + +\item{distribution}{the distribution used in the likelihood, see \link[complex]{clm}.} + +\item{h}{the forecast horizon.} + +\item{holdout}{if \code{TRUE}, the last \code{h} observations are held out.} + +\item{silent}{if \code{FALSE}, the progress is printed and the final model is plotted.} + +\item{...}{other parameters passed to \link[complex]{carima} and \link[complex]{clm}.} +} +\value{ +The selected model of the class \code{carima} (see \link[complex]{carima}), with +additional elements: +\itemize{ +\item \code{ICs} - data frame with the orders, the constant and the information criterion +(on the common sample) of all the models estimated in the process, +\item \code{ic} - the information criterion used, +\item \code{icObsDropped} - the number of first observations dropped for the criteria. +} +} +\description{ +Function selects the orders of complex ARIMA (see \link[complex]{carima}) based on an +information criterion. +} +\details{ +The selection follows the principles of \link[smooth]{auto.ssarima}: the AR and MA orders are +selected for each order of differencing, and the best model across all of them is chosen. +With \code{search="stepwise"}: +\enumerate{ +\item For each \eqn{d} from 0 to \code{orders$i}, the AR and MA orders are selected via a +neighbourhood search (similar to Hyndman & Khandakar, 2008): starting from cARIMA(0,d,0), +all models with p and q different by at most one are fitted, and the search moves to the best +of them, as long as the information criterion improves. If \code{constant=NULL}, the constant +is included for \eqn{d=0} and excluded otherwise; +\item If \code{constant=NULL}, the best model of each \eqn{d} is re-estimated with the opposite +option for the constant and the neighbourhood search continues from it. The best of all the +estimated models is returned. +} +With \code{search="full"}, all combinations of p, d, q and the constant are fitted +(except for cARIMA(p,0,q) without constant when \code{constant=NULL}). + +Models with different orders of differencing and AR use different numbers of initial values, +which \link[complex]{clm} constructs from the data. To make the information criteria +comparable, all models are estimated on the whole sample, but their criteria are calculated on +the common sample without the first \eqn{k} observations, where \eqn{k} is the sum of the +maximum orders. The log-likelihood on that sample is calculated at the estimated parameters, +with the scale re-estimated. + +Every model is estimated from the default starting values and from the parameters of the +closest model already estimated (with zeros for the new parameters), and the one with the +higher likelihood is kept. +} +\examples{ +set.seed(41) +e <- rcnorm(300, 0, sigma2=1, varsigma2=0.3+0.2i) +y <- cumsum(e + c(0, (0.5-0.2i)*e[-300])) +model <- auto.carima(y, orders=list(ar=2, i=1, ma=2)) +model +model$ICs + +} +\references{ +\itemize{ +\item Svetunkov, I., Boylan, J. E. (2019). State-space ARIMA for supply-chain forecasting. +International Journal of Production Research, 58(3), 818-827. +\item Hyndman, R. J., Khandakar, Y. (2008). Automatic time series forecasting: the forecast +package for R. Journal of Statistical Software, 27(3), 1-22. +} +} +\seealso{ +\code{\link[complex]{carima}}, \code{\link[smooth]{auto.ssarima}}, +\code{\link[smooth]{auto.msarima}} +} +\author{ +Ivan Svetunkov, \email{ivan@svetunkov.com} +} +\keyword{models} +\keyword{nonlinear} +\keyword{regression} +\keyword{ts} diff --git a/tests/testthat/test-auto-carima.R b/tests/testthat/test-auto-carima.R new file mode 100644 index 0000000..1fd1a8d --- /dev/null +++ b/tests/testthat/test-auto-carima.R @@ -0,0 +1,57 @@ +# auto.carima(): order selection for complex ARIMA + +test_that("auto.carima() finds cARIMA(1,1,1), the same as the full search", { + set.seed(41) + obs <- 400 + e <- rcnorm(obs, 0, sigma2=1, varsigma2=0.3+0.2i) + dy <- e + for(t in 2:obs){ + dy[t] <- (0.5+0.2i)*dy[t-1] + e[t] + (0.3-0.1i)*e[t-1] + } + y <- cumsum(dy) + model <- auto.carima(y, orders=list(ar=2, i=1, ma=2)) + expect_s3_class(model, "carima") + expect_equal(model$orders, list(ar=1, i=1, ma=1)) + expect_false(model$constant) + modelFull <- auto.carima(y, orders=list(ar=2, i=1, ma=2), search="full") + expect_equal(modelFull$orders, model$orders) + expect_equal(nrow(modelFull$ICs), 3*3*2 + 3*3) +}) + +test_that("auto.carima() finds the random walk with drift and the white noise", { + set.seed(7) + y <- cumsum(rcnorm(400, 0.3+0.1i, sigma2=1, varsigma2=0.2)) + model <- auto.carima(y, orders=list(ar=2, i=1, ma=2)) + expect_equal(model$orders, list(ar=0, i=1, ma=0)) + expect_true(model$constant) + y <- 5+2i + rcnorm(400, 0, sigma2=1, varsigma2=0.2) + model <- auto.carima(y, orders=list(ar=2, i=1, ma=2)) + expect_equal(model$orders, list(ar=0, i=0, ma=0)) + expect_true(model$constant) +}) + +test_that("the information criteria are calculated on the common sample", { + set.seed(41) + y <- cumsum(rcnorm(300, 0, sigma2=1, varsigma2=0.2)) + model <- auto.carima(y, orders=list(ar=1, i=1, ma=1), h=10, holdout=TRUE) + expect_equal(model$icObsDropped, 3) + expect_equal(nobs(model), 290) + expect_length(model$forecast, 10) + expect_true(!is.null(model$accuracy)) + # Without dropping observations the criterion would be the one of the model + modelFit <- carima(y[1:290], orders=c(model$orders$ar, model$orders$i, model$orders$ma), + constant=model$constant) + expect_equal(coef(modelFit), coef(model), tolerance=1e-4) +}) + +test_that("clm() penalises non-stationary and non-invertible parameters", { + set.seed(41) + y <- cumsum(rcnorm(300, 0, sigma2=1, varsigma2=0.2)) + model <- carima(y, orders=c(1,1,0), arma=list(ar=1.2)) + expect_lt(as.numeric(logLik(model)), -1e90) + model <- carima(y, orders=c(0,1,1), arma=list(ma=0.5-1i)) + expect_lt(as.numeric(logLik(model)), -1e90) + model <- carima(y, orders=c(1,1,1)) + expect_true(all(Mod(polyroot(c(1, -coef(model)[1])))>1)) + expect_true(all(Mod(polyroot(c(1, coef(model)[2])))>1)) +}) From 4f1fcafdbf10d26051ed2c4a085314bd3a3a95f0 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Sat, 26 Sep 2026 14:46:57 +0000 Subject: [PATCH 12/16] Fix clm() starting values for cARIMA(0,0,q) without intercept (1:0 indexing) Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NEWS | 1 + R/clm.R | 18 +++++++++--------- tests/testthat/test-arima-nointercept.R | 10 ++++++++++ 3 files changed, 20 insertions(+), 9 deletions(-) diff --git a/NEWS b/NEWS index 8b99125..0130122 100644 --- a/NEWS +++ b/NEWS @@ -15,6 +15,7 @@ Bugfixes: * clm() without intercept (e.g. y~-1) and with MA terms flipped the sign of the MA parameters (the likelihood was correct, the reported parameters and forecasts were not). * predict.clm() for cARIMA without intercept used the first parameter as an intercept (1:0 indexing). * clm() without intercept and with AR terms failed in the check of variability of differenced data. +* clm() with d=0, no intercept and no AR (e.g. cARIMA(0,0,q) with y~-1) selected the first column for the starting values (1:0 indexing) and failed. * clm() with nothing to estimate (e.g. cARIMA(0,1,0) without intercept) failed in the optimiser. * The documentation of the ellipsis in clm() was hidden by the one of sigma.clm(). diff --git a/R/clm.R b/R/clm.R index 6f0e93a..2b82d51 100644 --- a/R/clm.R +++ b/R/clm.R @@ -879,9 +879,9 @@ clm <- function(formula, data, subset, na.action, if(loss=="CLS"){ # If this is d=0 model if(iOrder==0){ - B <- as.vector(invert(t(matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE]) %*% - matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE]) %*% - t(matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE]) %*% y); + B <- as.vector(invert(t(matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE]) %*% + matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE]) %*% + t(matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE]) %*% y); } else{ matrixXregForDiffs <- iModelDesign(); @@ -903,9 +903,9 @@ clm <- function(formula, data, subset, na.action, else if(loss=="OLS"){ # If this is d=0 model if(iOrder==0){ - B <- as.vector(invert(t(Conj(matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE])) %*% - matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE]) %*% - t(Conj(matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE])) %*% y); + B <- as.vector(invert(t(Conj(matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE])) %*% + matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE]) %*% + t(Conj(matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE])) %*% y); } else{ matrixXregForDiffs <- iModelDesign(); @@ -932,9 +932,9 @@ clm <- function(formula, data, subset, na.action, if(is.null(B)){ # If this is d=0 model if(iOrder==0){ - B <- as.vector(invert(t(Conj(matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE])) %*% - matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE]) %*% - t(Conj(matrixXreg[,1:(nVariablesExo+arOrder), drop=FALSE])) %*% y); + B <- as.vector(invert(t(Conj(matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE])) %*% + matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE]) %*% + t(Conj(matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE])) %*% y); } else{ matrixXregForDiffs <- iModelDesign(); diff --git a/tests/testthat/test-arima-nointercept.R b/tests/testthat/test-arima-nointercept.R index 2226a54..681a659 100644 --- a/tests/testthat/test-arima-nointercept.R +++ b/tests/testthat/test-arima-nointercept.R @@ -44,3 +44,13 @@ test_that("cARIMA(0,1,0) without intercept has nothing to estimate and forecasts expect_equal(as.vector(predict(model, newdata=matrix(NA, 3, 1))$mean), rep(y[obs], 3)) } }) + +test_that("cARIMA(0,0,q) without intercept can be estimated", { + set.seed(4) + e <- rcnorm(500, 0, sigma2=1, varsigma2=0.3) + y <- e + c(0, (0.5-0.2i)*e[-500]) + for(loss in c("likelihood", "OLS", "CLS")){ + model <- suppressWarnings(clm(y~-1, data.frame(y=y), orders=c(0,0,1), loss=loss)) + expect_equal(unname(coef(model)), 0.5-0.2i, tolerance=0.1) + } +}) From b7ccaa434c0f27b87d931f338de7023d0208e810 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Sat, 26 Sep 2026 14:49:46 +0000 Subject: [PATCH 13/16] clm(): stop with an explanation for loss="CLS" with MA terms The CLS loss is complex-valued and cannot be minimised numerically, which the MA terms require; before, nloptr failed with an unclear error (also in 1.0.2). Fixes the test added in the previous commit. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NEWS | 1 + R/clm.R | 5 +++++ tests/testthat/test-arima-nointercept.R | 4 +++- 3 files changed, 9 insertions(+), 1 deletion(-) diff --git a/NEWS b/NEWS index 0130122..ac1f742 100644 --- a/NEWS +++ b/NEWS @@ -16,6 +16,7 @@ Bugfixes: * predict.clm() for cARIMA without intercept used the first parameter as an intercept (1:0 indexing). * clm() without intercept and with AR terms failed in the check of variability of differenced data. * clm() with d=0, no intercept and no AR (e.g. cARIMA(0,0,q) with y~-1) selected the first column for the starting values (1:0 indexing) and failed. +* clm() with loss="CLS" and MA terms failed with an unclear error from the optimiser. It now stops with an explanation (the CLS loss is complex-valued and cannot be minimised numerically). * clm() with nothing to estimate (e.g. cARIMA(0,1,0) without intercept) failed in the optimiser. * The documentation of the ellipsis in clm() was hidden by the one of sigma.clm(). diff --git a/R/clm.R b/R/clm.R index 2b82d51..77817c1 100644 --- a/R/clm.R +++ b/R/clm.R @@ -173,6 +173,11 @@ clm <- function(formula, data, subset, na.action, scaling <- match.arg(scaling); distribution <- match.arg(distribution); + if(loss=="CLS" && length(orders)==3 && orders[3]!=0){ + stop("loss=\"CLS\" cannot be used with MA terms: the parameters then need to be estimated ", + "numerically, while the CLS loss is complex-valued. Use loss=\"OLS\" or \"likelihood\".", + call.=FALSE); + } if(distribution=="dcgnorm" && loss!="likelihood"){ stop("distribution=\"dcgnorm\" can only be used with loss=\"likelihood\".", call.=FALSE); } diff --git a/tests/testthat/test-arima-nointercept.R b/tests/testthat/test-arima-nointercept.R index 681a659..f05ae9d 100644 --- a/tests/testthat/test-arima-nointercept.R +++ b/tests/testthat/test-arima-nointercept.R @@ -49,8 +49,10 @@ test_that("cARIMA(0,0,q) without intercept can be estimated", { set.seed(4) e <- rcnorm(500, 0, sigma2=1, varsigma2=0.3) y <- e + c(0, (0.5-0.2i)*e[-500]) - for(loss in c("likelihood", "OLS", "CLS")){ + for(loss in c("likelihood", "OLS")){ model <- suppressWarnings(clm(y~-1, data.frame(y=y), orders=c(0,0,1), loss=loss)) expect_equal(unname(coef(model)), 0.5-0.2i, tolerance=0.1) } + # CLS is complex-valued and cannot be used in the numeric optimisation needed for MA + expect_error(clm(y~-1, data.frame(y=y), orders=c(0,0,1), loss="CLS"), "cannot be used with MA") }) From 46f1ee0930bf425dda013266fda2b1a79951e762 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Sat, 26 Sep 2026 16:24:22 +0000 Subject: [PATCH 14/16] clm(fast=TRUE): Hannan-Rissanen starting values for pure cARIMA fast=TRUE (besides skipping the data checks) now starts the optimiser from Hannan-Rissanen estimates (long AR for the errors, then complex least squares on lags of the series and of the errors), shrunk to be admissible, and for dcgnorm from the shape maximising the likelihood of their residuals. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NEWS | 1 + R/clm.R | 39 ++++++++++- R/hannanRissanen.R | 93 +++++++++++++++++++++++++ man/clm.Rd | 6 +- tests/testthat/test-arima-nointercept.R | 18 +++++ 5 files changed, 154 insertions(+), 3 deletions(-) create mode 100644 R/hannanRissanen.R diff --git a/NEWS b/NEWS index ac1f742..979c093 100644 --- a/NEWS +++ b/NEWS @@ -5,6 +5,7 @@ Changes: * carima() - pure complex ARIMA (no regressors), a wrapper of clm() with an interface similar to smooth::msarima(): orders as c(p,d,q) or list(ar,i,ma), constant (FALSE by default, TRUE or a value), arma (fixed parameters), model (reuse of an estimated model), h, holdout (with accuracy), silent. Non-seasonal only. forecast() method for carima. * auto.carima() - automatic selection of the orders of complex ARIMA. Following auto.ssarima(), AR and MA orders are selected for each order of differencing via an IC-driven neighbourhood search (Hyndman & Khandakar, 2008), and the constant is tested at the end for the best model of each d. search="full" fits all the combinations. The information criteria are calculated on the common sample without the first ar+i+ma observations, so that models with different d are comparable. Each model is estimated from the default starting values and from the closest estimated model. * clm() now penalises non-stationary AR and non-invertible MA parameters in the loss (all roots of the complex AR/MA polynomials should lie outside the unit circle); the loss becomes 1E+100 divided by the smallest modulus of the roots. +* clm(fast=TRUE) for pure cARIMA (no explanatory variables, the intercept is allowed) now starts the optimiser from the Hannan-Rissanen estimates of AR and MA parameters (and, for dcgnorm, from the shape maximising the likelihood of their residuals). On wind data this reached much higher likelihoods than the default starting values (e.g. 489.8 vs 412.2 for cARIMA(1,1,2)). * dcgnorm() and rcgnorm() - density and random generation for the complex generalised normal distribution (complex generalised Gaussian), with scale, pseudo-scale and shape. * clm() has a new argument distribution. distribution="dcgnorm" uses the complex generalised normal likelihood, with the shape estimated (or provided via shape=..., with the starting value as the last element of B). The scale and the pseudo-scale are concentrated out: the circularity coefficient comes from the moments of the residuals and the scale is the ML estimate given it (approximate ML when shape is not 1). Standard errors are for the parameters given the shape. * AIC() and BIC() of clm() now count all real parameters of the model (both parts of the complex variable), as done in multivariate models. nparam() and the degrees of freedom are still per series. Before, the penalty was half of what it should be. Information criteria values change, and so can the orders selected by them. diff --git a/R/clm.R b/R/clm.R index 77817c1..0e4a6a9 100644 --- a/R/clm.R +++ b/R/clm.R @@ -58,7 +58,11 @@ #' is estimated. #' @param fast if \code{TRUE}, then the function won't check whether #' the data has variability and whether the regressors are correlated. Might -#' cause trouble, especially in cases of multicollinearity. +#' cause trouble, especially in cases of multicollinearity. In case of cARIMA without explanatory +#' variables (the intercept is allowed), \code{fast=TRUE} also makes the function use the +#' Hannan-Rissanen estimates (Hannan & Rissanen, 1982) of the AR and MA parameters as the starting +#' values of the optimiser (and, for \code{distribution="dcgnorm"}, the shape that maximises the +#' likelihood of their residuals), instead of the default ones. #' @param ... additional parameters to pass to distribution functions. This #' includes: #' \itemize{ @@ -879,6 +883,23 @@ clm <- function(formula, data, subset, na.action, # if(scaling!="none"){ # } + # Hannan-Rissanen starting values for pure cARIMA with fast=TRUE (NULL if not applicable or failed) + hrStart <- function(){ + if(!(fast && arimaModel && (arOrder>0 || maOrder>0) && (nVariablesExo-interceptIsNeeded)==0)){ + return(NULL); + } + u <- as.vector(if(iOrder>0) diff(y, differences=iOrder) else y); + estimates <- tryCatch(hannanRissanen(u, arOrder, maOrder, interceptIsNeeded), + error=function(e) NULL); + if(is.null(estimates)){ + return(NULL); + } + arValues <- hrAdmissible(estimates$ar, -1); + maValues <- hrAdmissible(estimates$ma, 1); + return(list(B=c(estimates$constant, arValues, maValues), + residuals=armaResiduals(u, estimates$constant, arValues, maValues))); + } + #### Estimate parameters of the model #### if(is.null(parameters)){ if(loss=="CLS"){ @@ -920,6 +941,10 @@ clm <- function(formula, data, subset, na.action, if(maOrderUsed){ # Add initial values for the maOrder B <- c(B, rep(0.1*(1+1i),maOrder)); + hrValues <- hrStart(); + if(!is.null(hrValues)){ + B <- hrValues$B; + } # Estimate the model res <- estimator(B, print_level); B <- res$B; @@ -934,7 +959,17 @@ clm <- function(formula, data, subset, na.action, # 1. paramExo, # 2. paramAR, # 3. paramMA. - if(is.null(B)){ + hrValues <- if(is.null(B)) hrStart() else NULL; + if(!is.null(hrValues)){ + B <- hrValues$B; + # The shape that maximises the likelihood of the Hannan-Rissanen residuals + if(shapeEstimate){ + hrResiduals <- hrValues$residuals[!is.na(hrValues$residuals)]; + shape <- exp(optimize(function(logShape) cgnormConcentrated(hrResiduals, exp(logShape))$logLik, + c(log(0.05), log(20)), maximum=TRUE)$maximum); + } + } + else if(is.null(B)){ # If this is d=0 model if(iOrder==0){ B <- as.vector(invert(t(Conj(matrixXreg[,seq_len(nVariablesExo+arOrder), drop=FALSE])) %*% diff --git a/R/hannanRissanen.R b/R/hannanRissanen.R new file mode 100644 index 0000000..a30d576 --- /dev/null +++ b/R/hannanRissanen.R @@ -0,0 +1,93 @@ +# Hannan-Rissanen estimation of complex ARMA(p,q) (Hannan & Rissanen, 1982), used for the +# starting values in clm(fast=TRUE) and for the screening in auto.carima(fast=TRUE). +# 1. A long AR(m) is estimated via complex least squares, and its residuals approximate the errors. +# 2. The series is regressed on its p lags and q lags of the approximate errors (and the constant), +# via complex least squares, which gives AR and MA parameters jointly. +# The estimates are consistent, but not efficient. + +# Complex least squares: minimises the sum of squared moduli of the residuals +complexLS <- function(X, target){ + XH <- t(Conj(X)); + return(as.vector(solve(XH %*% X, XH %*% target))); +} + +# Matrix of lags 1..k of the vector x (rows correspond to the elements of x, NA where unavailable) +lagMatrix <- function(x, k){ + n <- length(x); + X <- matrix(NA_complex_, n, k); + for(j in seq_len(k)){ + X[(j+1):n, j] <- x[seq_len(n-j)]; + } + return(X); +} + +# Residuals of ARMA(p,q) with the constant for the series u (all values available). +# The first p residuals are NA, the errors before the start are set to zero. +armaResiduals <- function(u, constant, ar, ma){ + n <- length(u); + p <- length(ar); + q <- length(ma); + errors <- rep(NA_complex_, n); + errorsLagged <- rep(0+0i, q); + constantValue <- if(is.null(constant)) 0 else constant; + for(t in (p+1):n){ + fitted <- constantValue; + if(p>0){ + fitted <- fitted + sum(ar * u[t-seq_len(p)]); + } + if(q>0){ + fitted <- fitted + sum(ma * errorsLagged); + } + errors[t] <- u[t] - fitted; + if(q>0){ + errorsLagged <- c(errors[t], errorsLagged[-q]); + } + } + return(errors); +} + +# Approximate errors from a long AR (stage 1). m is the order of the long AR. +hrErrors <- function(u, constant=FALSE, m=NULL){ + n <- length(u); + if(is.null(m)){ + m <- min(ceiling(10*log10(n)), floor((n-1)/4)); + } + X <- cbind(if(constant) rep(1+0i, n), lagMatrix(u, m)); + rows <- (m+1):n; + coefs <- complexLS(X[rows,,drop=FALSE], u[rows]); + errors <- rep(NA_complex_, n); + errors[rows] <- u[rows] - X[rows,,drop=FALSE] %*% coefs; + return(errors); +} + +# Stage 2: parameters of ARMA(p,q) with or without constant, given the approximate errors +hannanRissanen <- function(u, p, q, constant=FALSE, errors=NULL){ + n <- length(u); + if(p==0 && q==0){ + return(list(constant=if(constant) mean(u) else NULL, ar=NULL, ma=NULL)); + } + if(q>0 && is.null(errors)){ + errors <- hrErrors(u, constant); + } + X <- cbind(if(constant) rep(1+0i, n), lagMatrix(u, p), if(q>0) lagMatrix(errors, q)); + rows <- which(apply(!is.na(X), 1, all) & !is.na(u)); + coefs <- complexLS(X[rows,,drop=FALSE], u[rows]); + k <- as.integer(constant); + return(list(constant=if(constant) coefs[1] else NULL, + ar=coefs[k+seq_len(p)], + ma=coefs[k+p+seq_len(q)])); +} + +# Shrink AR / MA parameters towards zero until the roots of the polynomials are outside the unit circle +hrAdmissible <- function(parameters, sign){ + if(length(parameters)==0){ + return(parameters); + } + for(i in 1:100){ + if(all(Mod(polyroot(c(1, sign*parameters)))>1)){ + break; + } + parameters <- parameters*0.9; + } + return(parameters); +} diff --git a/man/clm.Rd b/man/clm.Rd index 41ea998..300303f 100644 --- a/man/clm.Rd +++ b/man/clm.Rd @@ -76,7 +76,11 @@ is estimated.} \item{fast}{if \code{TRUE}, then the function won't check whether the data has variability and whether the regressors are correlated. Might -cause trouble, especially in cases of multicollinearity.} +cause trouble, especially in cases of multicollinearity. In case of cARIMA without explanatory +variables (the intercept is allowed), \code{fast=TRUE} also makes the function use the +Hannan-Rissanen estimates (Hannan & Rissanen, 1982) of the AR and MA parameters as the starting +values of the optimiser (and, for \code{distribution="dcgnorm"}, the shape that maximises the +likelihood of their residuals), instead of the default ones.} \item{...}{additional parameters to pass to distribution functions. This includes: diff --git a/tests/testthat/test-arima-nointercept.R b/tests/testthat/test-arima-nointercept.R index f05ae9d..a5400f6 100644 --- a/tests/testthat/test-arima-nointercept.R +++ b/tests/testthat/test-arima-nointercept.R @@ -56,3 +56,21 @@ test_that("cARIMA(0,0,q) without intercept can be estimated", { # CLS is complex-valued and cannot be used in the numeric optimisation needed for MA expect_error(clm(y~-1, data.frame(y=y), orders=c(0,0,1), loss="CLS"), "cannot be used with MA") }) + +test_that("Hannan-Rissanen estimates and clm(fast=TRUE)", { + set.seed(41) + obs <- 4000 + e <- rcnorm(obs, 0, sigma2=1, varsigma2=0.3+0.2i) + u <- e + for(t in 2:obs){ + u[t] <- 1+0.5i + (0.5+0.2i)*u[t-1] + e[t] + (0.3-0.1i)*e[t-1] + } + estimates <- complex:::hannanRissanen(u, 1, 1, constant=TRUE) + expect_equal(c(estimates$constant, estimates$ar, estimates$ma), c(1+0.5i, 0.5+0.2i, 0.3-0.1i), + tolerance=0.1) + y <- cumsum(u) + model <- clm(y~1, data.frame(y=y), orders=c(1,1,1), fast=TRUE) + modelDefault <- clm(y~1, data.frame(y=y), orders=c(1,1,1)) + expect_gte(as.numeric(logLik(model)), as.numeric(logLik(modelDefault)) - 1e-3) + expect_equal(unname(coef(model)), c(1+0.5i, 0.5+0.2i, 0.3-0.1i), tolerance=0.1) +}) From f258cbb09bc140acc933f34e9d3edd7cbe07ae4f Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Sat, 26 Sep 2026 17:11:54 +0000 Subject: [PATCH 15/16] auto.carima(): Hannan-Rissanen starts by default; fast=TRUE screening Every candidate is estimated from the Hannan-Rissanen estimates (clm fast=TRUE) and from the closest estimated model. On wind data (log Z, 4 weeks) the selected cARIMA(1,1,1) with constant reached AICc -1040.93 vs -1034.15 with the default starts, with 55 fits instead of 61. fast=TRUE screens all models via Hannan-Rissanen and estimates via the likelihood the best one for each d, constant and MA order before the neighbourhood search: about three times faster, but approximate (on the same data it selected a model worse by 16 in AICc). Imports stats::optimize. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- NAMESPACE | 1 + NEWS | 1 + R/autoCarima.R | 193 ++++++++++++++++++++++++------ R/cgnorm.R | 2 +- man/auto.carima.Rd | 28 ++++- tests/testthat/test-auto-carima.R | 22 ++++ 6 files changed, 204 insertions(+), 43 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 88cb2a7..28d1c9c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -98,6 +98,7 @@ importFrom(stats, model.frame, model.matrix, nobs, + optimize, pacf, predict, qt, diff --git a/NEWS b/NEWS index 979c093..ca99438 100644 --- a/NEWS +++ b/NEWS @@ -4,6 +4,7 @@ complex v1.0.3 (Release data: 2026-09-25) Changes: * carima() - pure complex ARIMA (no regressors), a wrapper of clm() with an interface similar to smooth::msarima(): orders as c(p,d,q) or list(ar,i,ma), constant (FALSE by default, TRUE or a value), arma (fixed parameters), model (reuse of an estimated model), h, holdout (with accuracy), silent. Non-seasonal only. forecast() method for carima. * auto.carima() - automatic selection of the orders of complex ARIMA. Following auto.ssarima(), AR and MA orders are selected for each order of differencing via an IC-driven neighbourhood search (Hyndman & Khandakar, 2008), and the constant is tested at the end for the best model of each d. search="full" fits all the combinations. The information criteria are calculated on the common sample without the first ar+i+ma observations, so that models with different d are comparable. Each model is estimated from the default starting values and from the closest estimated model. +* auto.carima() starts every model from the Hannan-Rissanen estimates (via clm(fast=TRUE)) and from the closest estimated model. On wind data this reached a better optimum for the selected model (AICc lower by 6.8) with fewer fits. auto.carima(fast=TRUE) screens all the models via Hannan-Rissanen and estimates via the likelihood only the best one for each d, constant and MA order, before the neighbourhood search. It is about three times faster, but approximate. * clm() now penalises non-stationary AR and non-invertible MA parameters in the loss (all roots of the complex AR/MA polynomials should lie outside the unit circle); the loss becomes 1E+100 divided by the smallest modulus of the roots. * clm(fast=TRUE) for pure cARIMA (no explanatory variables, the intercept is allowed) now starts the optimiser from the Hannan-Rissanen estimates of AR and MA parameters (and, for dcgnorm, from the shape maximising the likelihood of their residuals). On wind data this reached much higher likelihoods than the default starting values (e.g. 489.8 vs 412.2 for cARIMA(1,1,2)). * dcgnorm() and rcgnorm() - density and random generation for the complex generalised normal distribution (complex generalised Gaussian), with scale, pseudo-scale and shape. diff --git a/R/autoCarima.R b/R/autoCarima.R index 3556ce7..e1bafc7 100644 --- a/R/autoCarima.R +++ b/R/autoCarima.R @@ -16,6 +16,17 @@ #' option for the constant and the neighbourhood search continues from it. The best of all the #' estimated models is returned. #' } +#' With \code{fast=TRUE}, all cARIMA(p,d,q) with and without constant are first estimated via +#' the Hannan-Rissanen method (Hannan & Rissanen, 1982): the differences are regressed on their +#' lags and on the lags of the errors approximated by a long AR, via complex least squares. Their +#' information criteria are approximate. Then the best of them for each \eqn{d}, constant and +#' MA order is estimated via the likelihood, and the neighbourhood search continues from the best +#' of these models (with its \eqn{d} and constant). With heavy tails, the least squares used in the +#' screening is not efficient and tends to underrate the models with MA terms, which is why the best +#' model for each MA order is estimated. This is still an approximation: on wind data with heavy +#' tails, \code{fast=TRUE} selected a model with an information criterion worse by 16 than +#' the one found with \code{fast=FALSE}, while being about three times faster. +#' #' With \code{search="full"}, all combinations of p, d, q and the constant are fitted #' (except for cARIMA(p,0,q) without constant when \code{constant=NULL}). #' @@ -26,9 +37,9 @@ #' maximum orders. The log-likelihood on that sample is calculated at the estimated parameters, #' with the scale re-estimated. #' -#' Every model is estimated from the default starting values and from the parameters of the -#' closest model already estimated (with zeros for the new parameters), and the one with the -#' higher likelihood is kept. +#' Every model is estimated from the Hannan-Rissanen estimates of its parameters (see +#' \code{fast} in \link[complex]{clm}) and from the parameters of the closest model already +#' estimated (with zeros for the new parameters), and the one with the higher likelihood is kept. #' #' @template author #' @template keywords @@ -44,6 +55,8 @@ #' @param h the forecast horizon. #' @param holdout if \code{TRUE}, the last \code{h} observations are held out. #' @param silent if \code{FALSE}, the progress is printed and the final model is plotted. +#' @param fast if \code{TRUE}, the models are screened via the Hannan-Rissanen method before the +#' search (see details). Ignored with \code{search="full"}. #' @param ... other parameters passed to \link[complex]{carima} and \link[complex]{clm}. #' #' @return The selected model of the class \code{carima} (see \link[complex]{carima}), with @@ -52,10 +65,14 @@ #' \item \code{ICs} - data frame with the orders, the constant and the information criterion #' (on the common sample) of all the models estimated in the process, #' \item \code{ic} - the information criterion used, -#' \item \code{icObsDropped} - the number of first observations dropped for the criteria. +#' \item \code{icObsDropped} - the number of first observations dropped for the criteria, +#' \item \code{ICsScreening} - with \code{fast=TRUE}, the approximate criteria of all the models +#' estimated via Hannan-Rissanen. #' } #' #' @references \itemize{ +#' \item Hannan, E. J., Rissanen, J. (1982). Recursive estimation of mixed autoregressive-moving +#' average order. Biometrika, 69(1), 81-94. #' \item Svetunkov, I., Boylan, J. E. (2019). State-space ARIMA for supply-chain forecasting. #' International Journal of Production Research, 58(3), 818-827. #' \item Hyndman, R. J., Khandakar, Y. (2008). Automatic time series forecasting: the forecast @@ -78,7 +95,7 @@ auto.carima <- function(y, orders=list(ar=3, i=2, ma=3), constant=NULL, ic=c("AICc","AIC","BIC","BICc"), search=c("stepwise","full"), distribution=c("dcnorm","dcgnorm"), - h=0, holdout=FALSE, silent=TRUE, ...){ + h=0, holdout=FALSE, silent=TRUE, fast=FALSE, ...){ cl <- match.call(); ic <- match.arg(ic); search <- match.arg(search); @@ -86,6 +103,10 @@ auto.carima <- function(y, orders=list(ar=3, i=2, ma=3), constant=NULL, ellipsis <- list(...); # The starting values are set by the function ellipsis$B <- NULL; + if(fast && search=="full"){ + warning("fast=TRUE is only used with search=\"stepwise\". Ignoring it.", call.=FALSE); + fast <- FALSE; + } orders <- carimaOrders(orders); arMax <- orders$ar; @@ -153,10 +174,12 @@ auto.carima <- function(y, orders=list(ar=3, i=2, ma=3), constant=NULL, if(!is.null(estimated[[key]])){ return(estimated[[key]]); } + # Without B, clm(fast=TRUE) starts from the Hannan-Rissanen estimates fit <- function(B=NULL){ carimaArgs <- c(list(y=y, orders=c(p, d, q), constant=constant, loss="likelihood", distribution=distribution), - if(!is.null(B)) list(B=B), ellipsis); + if(!is.null(B)) list(B=B) else list(fast=TRUE), + ellipsis); tryCatch(do.call(carima, carimaArgs), error=function(e) NULL); } fits <- c(list(fit()), @@ -180,7 +203,128 @@ auto.carima <- function(y, orders=list(ar=3, i=2, ma=3), constant=NULL, constantOptions <- if(is.null(constant)) c(FALSE, TRUE) else constant; - if(search=="full"){ + #### Screening via Hannan-Rissanen (fast=TRUE) #### + # For each d and constant, all cARIMA(p,d,q) are estimated via Hannan-Rissanen (complex least + # squares on the differences, with the errors approximated by a long AR). Their information + # criteria (from the residuals of the recursion, on the common sample) are approximate. + obsInSample <- length(y); + # Information criterion from the residuals (aligned with y) on the common sample, and the + # shape of dcgnorm that maximises the likelihood of these residuals + screenIC <- function(errors, nComplex){ + errors <- errors[(obsDropped+1):obsInSample]; + obsTrimmed <- length(errors); + shapeEstimated <- distribution=="dcgnorm"; + shapeValue <- NULL; + if(shapeEstimated){ + shapeFit <- optimize(function(logShape) cgnormConcentrated(errors, exp(logShape))$logLik, + c(log(0.05), log(20)), maximum=TRUE); + llikelihood <- shapeFit$objective; + shapeValue <- exp(shapeFit$maximum); + nParam <- nComplex + 3/2 + 1/2; + } + else{ + errorsMatrix <- complex2vec(errors); + llikelihood <- -obsTrimmed*(log(2*pi) + 1 + + 0.5*log(det(t(errorsMatrix) %*% errorsMatrix / obsTrimmed))); + nParam <- nComplex + 3/2; + } + return(list(IC=clmIC(llikelihood, nParam, obsTrimmed, shapeEstimated, ic), shape=shapeValue)); + } + hrScreen <- function(){ + screenRows <- list(); + for(d in 0:iMax){ + for(constantValue in constantOptions){ + if(d==0 && !constantValue && is.null(constant)){ + next; + } + u <- as.vector(if(d>0) diff(y, differences=d) else y); + errors <- if(maMax>0) tryCatch(hrErrors(u, constantValue), error=function(e) NULL) else NULL; + for(p in 0:arMax){ + for(q in 0:maMax){ + key <- modelKey(p, d, q, constantValue); + ICValue <- Inf; + estimates <- NULL; + if(q==0 || !is.null(errors)){ + estimates <- tryCatch(hannanRissanen(u, p, q, constantValue, errors), + error=function(e) NULL); + } + if(!is.null(estimates) && + (length(estimates$ar)==0 || all(Mod(polyroot(c(1, -estimates$ar)))>1)) && + (length(estimates$ma)==0 || all(Mod(polyroot(c(1, estimates$ma)))>1))){ + residualsHR <- c(rep(NA_complex_, d), + armaResiduals(u, estimates$constant, estimates$ar, estimates$ma)); + screenValues <- screenIC(residualsHR, as.integer(constantValue) + p + q); + ICValue <- screenValues$IC; + } + screenRows[[key]] <- data.frame(p=p, d=d, q=q, constant=constantValue, + IC=ICValue, key=key); + if(!silent){ + cat("Hannan-Rissanen cARIMA(", p, ",", d, ",", q, ")", + ifelse(constantValue, " with constant", ""), ": ", round(ICValue, 3), + "\n", sep=""); + } + } + } + } + } + screenTable <- do.call(rbind, screenRows); + rownames(screenTable) <- NULL; + return(screenTable); + } + + # Neighbourhood search: from the given model, move to the best of the models with p and q + # different by at most one, as long as the information criterion improves + neighbourhoodSearch <- function(best, d, constantValue){ + pBest <- best$model$orders$ar; + qBest <- best$model$orders$ma; + repeat{ + neighbours <- expand.grid(p=pBest+(-1:1), q=qBest+(-1:1)); + neighbours <- neighbours[neighbours$p>=0 & neighbours$p<=arMax & + neighbours$q>=0 & neighbours$q<=maMax & + !(neighbours$p==pBest & neighbours$q==qBest),,drop=FALSE]; + if(nrow(neighbours)==0){ + break; + } + results <- lapply(seq_len(nrow(neighbours)), function(i){ + fitCandidate(neighbours$p[i], d, neighbours$q[i], constantValue, + if(!is.null(best$model)) list(best$model) else list()); + }); + neighbourICs <- sapply(results, `[[`, "IC"); + if(min(neighbourICs) < best$IC){ + best <- results[[which.min(neighbourICs)]]; + pBest <- neighbours$p[which.min(neighbourICs)]; + qBest <- neighbours$q[which.min(neighbourICs)]; + } + else{ + break; + } + } + return(best); + } + + if(fast){ + screenTable <- hrScreen(); + # The best screened model for each d, constant and q is estimated via the likelihood. + # Taking the best for each q protects the models with MA terms, which the screening can + # underrate (e.g. with heavy tails, where least squares is not efficient) + screenTable <- screenTable[order(screenTable$IC),,drop=FALSE]; + bestScreened <- screenTable[is.finite(screenTable$IC) & + !duplicated(screenTable[,c("d","constant","q")]),,drop=FALSE]; + if(!silent){ + cat("Estimating", nrow(bestScreened), "models via the likelihood...\n"); + } + for(i in seq_len(nrow(bestScreened))){ + fitCandidate(bestScreened$p[i], bestScreened$d[i], bestScreened$q[i], bestScreened$constant[i]); + } + # The neighbourhood search with the likelihood continues from the best of them + ICs <- sapply(estimated, `[[`, "IC"); + best <- estimated[[which.min(ICs)]]; + if(!silent){ + cat("Neighbourhood search from the best model...\n"); + } + neighbourhoodSearch(best, best$model$orders$i, best$model$constant); + } + else if(search=="full"){ for(d in 0:iMax){ for(constantValue in constantOptions){ if(d==0 && !constantValue && is.null(constant)){ @@ -207,36 +351,6 @@ auto.carima <- function(y, orders=list(ar=3, i=2, ma=3), constant=NULL, } } else{ - # Neighbourhood search: from the given model, move to the best of the models with p and q - # different by at most one, as long as the information criterion improves - neighbourhoodSearch <- function(best, d, constantValue){ - pBest <- best$model$orders$ar; - qBest <- best$model$orders$ma; - repeat{ - neighbours <- expand.grid(p=pBest+(-1:1), q=qBest+(-1:1)); - neighbours <- neighbours[neighbours$p>=0 & neighbours$p<=arMax & - neighbours$q>=0 & neighbours$q<=maMax & - !(neighbours$p==pBest & neighbours$q==qBest),,drop=FALSE]; - if(nrow(neighbours)==0){ - break; - } - results <- lapply(seq_len(nrow(neighbours)), function(i){ - fitCandidate(neighbours$p[i], d, neighbours$q[i], constantValue, - if(!is.null(best$model)) list(best$model) else list()); - }); - neighbourICs <- sapply(results, `[[`, "IC"); - if(min(neighbourICs) < best$IC){ - best <- results[[which.min(neighbourICs)]]; - pBest <- neighbours$p[which.min(neighbourICs)]; - qBest <- neighbours$q[which.min(neighbourICs)]; - } - else{ - break; - } - } - return(best); - } - # Search for p and q within each d, starting from cARIMA(0,d,0). The constant is included # for d=0 and excluded for d>0 (unless it is provided), and is tested at the end. bestPerDiff <- vector("list", iMax+1); @@ -289,6 +403,11 @@ auto.carima <- function(y, orders=list(ar=3, i=2, ma=3), constant=NULL, rownames(ICsTable) <- NULL; bestModel$ICs <- ICsTable; + if(fast){ + screenTable <- screenTable[, c("p","d","q","constant","IC")]; + names(screenTable)[5] <- ic; + bestModel$ICsScreening <- screenTable; + } bestModel$ic <- ic; bestModel$icObsDropped <- obsDropped; bestModel$autoCall <- cl; diff --git a/R/cgnorm.R b/R/cgnorm.R index e7a66bb..18701b5 100644 --- a/R/cgnorm.R +++ b/R/cgnorm.R @@ -73,7 +73,7 @@ dcgnorm <- function(q, mu=0, scale=1, pseudoscale=0, shape=1, log=FALSE){ #' @rdname cgnormal #' @export rcgnorm #' @aliases rcgnorm -#' @importFrom stats rgamma runif +#' @importFrom stats rgamma runif optimize rcgnorm <- function(n=1, mu=0, scale=1, pseudoscale=0, shape=1){ cgnormCheck(scale, pseudoscale, shape); # Covariance-type matrix of the real and imaginary parts (as in dcnorm) diff --git a/man/auto.carima.Rd b/man/auto.carima.Rd index cf0cfae..dad889f 100644 --- a/man/auto.carima.Rd +++ b/man/auto.carima.Rd @@ -7,7 +7,7 @@ auto.carima(y, orders = list(ar = 3, i = 2, ma = 3), constant = NULL, ic = c("AICc", "AIC", "BIC", "BICc"), search = c("stepwise", "full"), distribution = c("dcnorm", "dcgnorm"), h = 0, holdout = FALSE, - silent = TRUE, ...) + silent = TRUE, fast = FALSE, ...) } \arguments{ \item{y}{vector or ts object with the complex-valued time series.} @@ -30,6 +30,9 @@ combinations of orders.} \item{silent}{if \code{FALSE}, the progress is printed and the final model is plotted.} +\item{fast}{if \code{TRUE}, the models are screened via the Hannan-Rissanen method before the +search (see details). Ignored with \code{search="full"}.} + \item{...}{other parameters passed to \link[complex]{carima} and \link[complex]{clm}.} } \value{ @@ -39,7 +42,9 @@ additional elements: \item \code{ICs} - data frame with the orders, the constant and the information criterion (on the common sample) of all the models estimated in the process, \item \code{ic} - the information criterion used, -\item \code{icObsDropped} - the number of first observations dropped for the criteria. +\item \code{icObsDropped} - the number of first observations dropped for the criteria, +\item \code{ICsScreening} - with \code{fast=TRUE}, the approximate criteria of all the models +estimated via Hannan-Rissanen. } } \description{ @@ -60,6 +65,17 @@ is included for \eqn{d=0} and excluded otherwise; option for the constant and the neighbourhood search continues from it. The best of all the estimated models is returned. } +With \code{fast=TRUE}, all cARIMA(p,d,q) with and without constant are first estimated via +the Hannan-Rissanen method (Hannan & Rissanen, 1982): the differences are regressed on their +lags and on the lags of the errors approximated by a long AR, via complex least squares. Their +information criteria are approximate. Then the best of them for each \eqn{d}, constant and +MA order is estimated via the likelihood, and the neighbourhood search continues from the best +of these models (with its \eqn{d} and constant). With heavy tails, the least squares used in the +screening is not efficient and tends to underrate the models with MA terms, which is why the best +model for each MA order is estimated. This is still an approximation: on wind data with heavy +tails, \code{fast=TRUE} selected a model with an information criterion worse by 16 than +the one found with \code{fast=FALSE}, while being about three times faster. + With \code{search="full"}, all combinations of p, d, q and the constant are fitted (except for cARIMA(p,0,q) without constant when \code{constant=NULL}). @@ -70,9 +86,9 @@ the common sample without the first \eqn{k} observations, where \eqn{k} is the s maximum orders. The log-likelihood on that sample is calculated at the estimated parameters, with the scale re-estimated. -Every model is estimated from the default starting values and from the parameters of the -closest model already estimated (with zeros for the new parameters), and the one with the -higher likelihood is kept. +Every model is estimated from the Hannan-Rissanen estimates of its parameters (see +\code{fast} in \link[complex]{clm}) and from the parameters of the closest model already +estimated (with zeros for the new parameters), and the one with the higher likelihood is kept. } \examples{ set.seed(41) @@ -85,6 +101,8 @@ model$ICs } \references{ \itemize{ +\item Hannan, E. J., Rissanen, J. (1982). Recursive estimation of mixed autoregressive-moving +average order. Biometrika, 69(1), 81-94. \item Svetunkov, I., Boylan, J. E. (2019). State-space ARIMA for supply-chain forecasting. International Journal of Production Research, 58(3), 818-827. \item Hyndman, R. J., Khandakar, Y. (2008). Automatic time series forecasting: the forecast diff --git a/tests/testthat/test-auto-carima.R b/tests/testthat/test-auto-carima.R index 1fd1a8d..5787734 100644 --- a/tests/testthat/test-auto-carima.R +++ b/tests/testthat/test-auto-carima.R @@ -55,3 +55,25 @@ test_that("clm() penalises non-stationary and non-invertible parameters", { expect_true(all(Mod(polyroot(c(1, -coef(model)[1])))>1)) expect_true(all(Mod(polyroot(c(1, coef(model)[2])))>1)) }) + +test_that("auto.carima(fast=TRUE) finds the same models via Hannan-Rissanen screening", { + set.seed(41) + obs <- 400 + e <- rcnorm(obs, 0, sigma2=1, varsigma2=0.3+0.2i) + dy <- e + for(t in 2:obs){ + dy[t] <- (0.5+0.2i)*dy[t-1] + e[t] + (0.3-0.1i)*e[t-1] + } + y <- cumsum(dy) + model <- auto.carima(y, orders=list(ar=2, i=1, ma=2), fast=TRUE) + expect_equal(model$orders, list(ar=1, i=1, ma=1)) + expect_false(model$constant) + expect_true(!is.null(model$ICsScreening)) + expect_equal(nrow(model$ICsScreening), 3*3*3) + set.seed(7) + y <- cumsum(rcnorm(400, 0.3+0.1i, sigma2=1, varsigma2=0.2)) + model <- auto.carima(y, orders=list(ar=2, i=1, ma=2), fast=TRUE) + expect_equal(model$orders, list(ar=0, i=1, ma=0)) + expect_true(model$constant) + expect_warning(auto.carima(y, orders=list(ar=1, i=1, ma=1), search="full", fast=TRUE), "fast") +}) From c093f8b608c925acae736b9219153eb98ad30715 Mon Sep 17 00:00:00 2001 From: Ivan Svetunkov Date: Sat, 26 Sep 2026 18:36:13 +0000 Subject: [PATCH 16/16] Add CLAUDE.md with the package layout, workflow and internals Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01V6iwpznpavEsZmK26Xf51m --- .Rbuildignore | 1 + CLAUDE.md | 92 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 93 insertions(+) create mode 100644 CLAUDE.md diff --git a/.Rbuildignore b/.Rbuildignore index 9048bdd..d46ba46 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -10,3 +10,4 @@ man-roxygen cran-comments.md ^CRAN-SUBMISSION$ ^dev$ +^CLAUDE\.md$ diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..3394532 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,92 @@ +# complex (R package) + +Complex-valued econometrics and forecasting: complex linear models, complex ARIMA (cARIMA), complex +distributions, complex correlations. Companion to Svetunkov & Svetunkov (2024), *Complex-Valued +Econometrics with Examples in R* (Springer); Chapter 6 covers the dynamic models. +Author/maintainer: Ivan Svetunkov. Depends on greybox and legion (which loads smooth). + +## Repository + +- `git@github.com:openforecast-org/complex.git` (moved from config-i1/complex). Always use ssh. + There is no ssh in the Claude container, so the user does push/fetch. +- Branch `dcgnorm`: development of v1.0.3 (not yet merged to `master`). Commit as + `Ivan Svetunkov ` (`git -c user.name=... -c user.email=...`). +- Used by the wind forecasting paper project: `~/R/Projects/Experiments/cARIMA` + (git@github.com:config-i1/complex-wind.git), which has its own CLAUDE.md. + +## Layout + +| file | content | +|---|---| +| `R/clm.R` | `clm()` (complex linear model and cARIMA(X)) and its methods: `logLik`, `nparam`, `AICc/BICc`, `clmIC()`, `vcov`, `confint`, `summary`, `predict`, `plot` | +| `R/carima.R` | `carima()` pure cARIMA wrapper of `clm()` (msarima-like interface), `forecast.carima`, `carimaFinalise()` | +| `R/autoCarima.R` | `auto.carima()` order selection (stepwise / full / fast) | +| `R/hannanRissanen.R` | Hannan-Rissanen estimation (internal): starting values and screening | +| `R/cnorm.R`, `R/cgnorm.R` | complex normal; complex generalised normal (`dcgnorm`, `rcgnorm`, internal `cgnormConcentrated`) | +| `R/cacf.R`, `R/cvar.R` | complex ACF/PACF, direct/conjugate variances and correlations | +| `R/cTransformations.R` | `cscale`, `cdescale`, `clog`, `cexp` | +| `src/complexCode.cpp` | Rcpp: `invert`, `polyprodcomplex` | +| `tests/testthat/` | tests (none existed before v1.0.3) | +| `dev/` | plans and prototypes (in `.Rbuildignore`); `dev/PLAN-dcgnorm.md` has the design decisions | + +## Workflow + +```bash +cd ~/R/Projects/Packages/complex +Rscript -e 'roxygen2::roxygenise()' # roxygen2 8.1.0; never edit NAMESPACE/man by hand +R CMD INSTALL . # installs into the user library +Rscript -e 'library(complex); testthat::test_dir("tests/testthat", reporter="summary")' +R CMD build --no-build-vignettes . && R CMD check --no-manual --ignore-vignettes complex_*.tar.gz +``` +- Build and check in a scratch directory, not in the package folder. +- Expected check result: `Status: 1 WARNING`, the "OS reports request to set locale to en_US.UTF-8" + warning of this machine. Anything else must be fixed. +- Add a NEWS entry (Changes / Bugfixes) for every user-visible change. +- `src/*.o` and `src/complex.so` are tracked in git (historical). Do not commit changes to them. + +## Code style + +Match the existing code: 4-space indentation, `;` at the end of statements, `if(...){ ... }` with +`else{` on the next line after `}`, camelCase names, `<-`, comments describing why. Use +`seq_len(n)` rather than `1:n`: several bugs here came from `1:0` selecting the first element. + +## clm() internals worth knowing + +- Parameters `B` are complex (intercept, regressors, AR, MA). The optimiser (nloptr, SBPLX by default) + works on `c(Re(B), Im(B))`; for `distribution="dcgnorm"` with estimated shape, log(shape) is + appended at the end. A starting shape can be passed as the last element of `B`. +- The ARI part is stored via the polynomial `(1-phi(B))(1-B)^d`; `fitter()` builds the full + parameter vector. Initial lags are extrapolated (`xregExpander(gaps="auto")`), initial MA errors + are zero; all observations enter the likelihood. +- Likelihood `dcnorm`: concentrated, `-T(log 2pi + 1 + 0.5 log det Sigma)`. `dcgnorm`: circularity + from moments, scale by exact ML given it (approximate ML for shape != 1). +- The loss returns `1E+100/min|root|` if the complex AR or MA polynomial has a root inside the unit + circle (differencing is not checked). +- `nparam()` counts per series (complex units: covariance 3/2, shape 1/2). `logLik` df and the ICs + use all real parameters (2*nparam); `clmIC()` has the formulas (AICc/BICc: Bedrick & Tsai with + two series, approximate for the restricted complex model). +- `fast=TRUE`: skips data checks and, for pure cARIMA, starts from Hannan-Rissanen estimates. +- `vcov.clm` re-estimates the model from `object$call` with `parameters` fixed and `FI=TRUE`; the + MA columns of `object$data` then act as regressors, so the parameter vector does not match the + orders there. Keep `object$call` a clm call (carima stores its own call in `carimaCall`). +- `loss="CLS"` cannot be used with MA (complex-valued loss). + +## auto.carima() + +- stepwise (default): for each d, neighbourhood search over p, q (IC only); constant included for + d=0, excluded for d>0, then flipped for the best model of each d and the search continued. +- Every candidate starts from Hannan-Rissanen and from the closest estimated model. +- ICs are computed on the common sample without the first `ar+i+ma` observations (models are + estimated on the full sample). +- `fast=TRUE`: Hannan-Rissanen screening of all models, likelihood refit of the best per + (d, constant, q), then neighbourhood search. About 3x faster, but approximate (least squares + underrates MA models with heavy tails). +- `search="full"`: all combinations. + +## Open items + +- Phase 3/4/5 of `dev/PLAN-dcgnorm.md`: joint Hessian with the shape, simulated prediction + intervals, dcgnorm diagnostics. +- Seasonal cARIMA is not supported. +- A complex portmanteau / partial autocorrelation statistic combining direct and conjugate ACFs + (discussed, not implemented).