diff --git a/.Rbuildignore b/.Rbuildignore index 9fb1e4e..d46ba46 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -9,3 +9,5 @@ man-roxygen ^revdep$ 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). diff --git a/DESCRIPTION b/DESCRIPTION index 67e1335..42de076 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,13 @@ 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 -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 @@ -27,6 +28,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 7d94243..28d1c9c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,7 +1,10 @@ # Generated by roxygen2: do not edit by hand +S3method(AICc,clm) +S3method(BICc,clm) S3method(actuals,clm) S3method(confint,clm) +S3method(forecast,carima) S3method(logLik,clm) S3method(nobs,clm) S3method(nparam,clm) @@ -10,11 +13,15 @@ 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(auto.carima) export(cacf) +export(carima) export(ccor) export(ccov) export(ccov2cor) @@ -29,70 +36,89 @@ 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) -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, + forecast, + 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, + optimize, + pacf, + predict, + qt, + resid, + residuals, + rgamma, + rnorm, + runif, + sd, + setNames, + sigma, + start, + terms, + time, + ts, + update.formula, + vcov +) +importFrom(utils, + packageVersion, + tail +) useDynLib(complex) diff --git a/NEWS b/NEWS index 4eaee1d..ca99438 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,28 @@ +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. +* 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 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(). + + complex v1.0.2 (Release data: 2026-02-03) ============== diff --git a/R/autoCarima.R b/R/autoCarima.R new file mode 100644 index 0000000..e1bafc7 --- /dev/null +++ b/R/autoCarima.R @@ -0,0 +1,417 @@ +#' 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{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}). +#' +#' 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 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 +#' +#' @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 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 +#' 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, +#' \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 +#' 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, fast=FALSE, ...){ + 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; + 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; + 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]]); + } + # 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) else list(fast=TRUE), + 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; + + #### 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)){ + 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{ + # 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; + 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; + bestModel <- carimaFinalise(bestModel, h, holdout, holdoutData, silent); + + return(bestModel); +} diff --git a/R/carima.R b/R/carima.R new file mode 100644 index 0000000..46e4736 --- /dev/null +++ b/R/carima.R @@ -0,0 +1,254 @@ +#' 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)); + + 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){ + 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/R/cgnorm.R b/R/cgnorm.R new file mode 100644 index 0000000..18701b5 --- /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 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) + 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 95f2a87..0e4a6a9 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. @@ -48,13 +58,19 @@ #' 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{ #' \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: @@ -62,6 +78,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 @@ -136,6 +154,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 +176,15 @@ 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); + } #### Functions used in the estimation #### ifelseFast <- function(condition, yes, no){ @@ -197,7 +225,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 @@ -218,7 +246,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); } } @@ -257,9 +286,39 @@ 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)]; + } + + # 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"){ + 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 +380,16 @@ 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; + } + + # 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, @@ -329,7 +398,11 @@ clm <- function(formula, data, subset, na.action, # lb=BLower, ub=BUpper, loss=loss, y=y, matrixXreg=matrixXreg); BReal[] <- res$solution; - B[] <- complex(real=BReal[1:(nVariables/2)],imaginary=BReal[(nVariables/2+1):nVariables]); + shapeValue <- shape; + if(shapeEstimate){ + shapeValue <- exp(BReal[nVariables+1]); + } + B[] <- complex(real=BReal[seq_len(nVariables/2)],imaginary=BReal[nVariables/2+seq_len(nVariables/2)]); nVariables <- length(B); CFValue <- res$objective; @@ -337,13 +410,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; @@ -766,9 +858,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)){ @@ -776,7 +873,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) @@ -786,14 +883,31 @@ 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"){ # 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(); @@ -815,9 +929,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(); @@ -827,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; @@ -841,12 +959,22 @@ 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[,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(); @@ -858,9 +986,16 @@ 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; + shape <- res$shape; } } # If the parameters are provided @@ -895,6 +1030,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 +1063,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 +1101,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, @@ -994,7 +1137,50 @@ 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). The correction is derived for unrestricted +# multivariate regression, so it is approximate for clm. +#' @export +AICc.clm <- function(object, ...){ + return(clmIC(logLik(object), nparam(object), nobs(object), + isTRUE(object$other$shapeEstimated), "AICc")); +} + +#' @export +BICc.clm <- function(object, ...){ + 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; + 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){ + return(Inf); + } + penalty <- obs*nParamAll/(obs - (nParamPerSeries + nSeries + 1)); + if(ic=="AICc"){ + return(-2*llikelihood + 2*penalty); + } + return(-2*llikelihood + log(obs)*penalty); } #' @rdname clm @@ -1005,7 +1191,6 @@ logLik.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, ...){ @@ -1173,6 +1358,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 @@ -1294,6 +1484,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]]; @@ -1325,6 +1516,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)), @@ -1407,7 +1603,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); @@ -1501,10 +1698,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/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/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") ``` diff --git a/dev/PLAN-dcgnorm.md b/dev/PLAN-dcgnorm.md new file mode 100644 index 0000000..a1d3db8 --- /dev/null +++ b/dev/PLAN-dcgnorm.md @@ -0,0 +1,162 @@ +# 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 (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) + +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. + +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 ρ. + +## 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. +Commit per phase. + +### Phase 1 — distribution functions (`R/cgnorm.R`, new) +- `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, 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")` (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 + 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`. +- `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. 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: 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); 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`. 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") diff --git a/man/auto.carima.Rd b/man/auto.carima.Rd new file mode 100644 index 0000000..dad889f --- /dev/null +++ b/man/auto.carima.Rd @@ -0,0 +1,122 @@ +% 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, fast = FALSE, ...) +} +\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{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{ +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, +\item \code{ICsScreening} - with \code{fast=TRUE}, the approximate criteria of all the models +estimated via Hannan-Rissanen. +} +} +\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{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}). + +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 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) +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 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 +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/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/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..300303f 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. @@ -64,9 +76,45 @@ 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: +\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; +} -\item{...}{Other parameters passed to internal functions.} +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.} 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-arima-nointercept.R b/tests/testthat/test-arima-nointercept.R new file mode 100644 index 0000000..a5400f6 --- /dev/null +++ b/tests/testthat/test-arima-nointercept.R @@ -0,0 +1,76 @@ +# 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)) + } +}) + +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")){ + 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") +}) + +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) +}) diff --git a/tests/testthat/test-auto-carima.R b/tests/testthat/test-auto-carima.R new file mode 100644 index 0000000..5787734 --- /dev/null +++ b/tests/testthat/test-auto-carima.R @@ -0,0 +1,79 @@ +# 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)) +}) + +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") +}) 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) +}) 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")) +}) 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))) +})