R/fit.R
In GeoBosh/sarima: Simulation and Prediction with Seasonal ARIMA Models
Defines functions
sarimat
sarima
Documented in
sarima
## Do not edit this file manually.
## It has been automatically generated from *.org sources.
sarima <- function(model, data = NULL, ss.method = "sarima", use.symmetry = FALSE,
SSinit = "Rossignol2011"){
# TODO: for now lik.method is simply "ss.method"
# Need to think over the user-facing arguments for methods.
## N
lik.method <- ss.method
if(!inherits(model, "formula"))
stop("unsupported class ", class(model), " of argument 'model'")
toPolyCoef <- function(item){
switch(item$name,
ar = { phi <<- c( phi , list(item$coef) ) },
ma = { theta <<- c( theta , list(item$coef) ) },
sar = { phi <<- c( phi , list(item$coef) ) },
sma = { theta <<- c( theta , list(item$coef) ) },
i = { delta <<- c( delta , list(item$coef) ) },
si = { delta <<- c( delta , list(item$coef) ) },
s = { delta <<- c( delta , list(item$coef) ) },
ss = { delta <<- c( delta , list(item$coef) ) },
# here item$coef is itself a list!
su = { delta <<- c( delta , item$coef ) },
u = { delta <<- c( delta , item$coef ) },
# unit root for estimation:
uar = { udelta <<- c( udelta, list(item$coef) ) },
stop("not supported yet!")
)
}
Fo <- Formula(model)
yname <- all.vars(formula(Fo, rhs = FALSE)) # name of the response variable
if(length(yname) != 1)
stop("currently there should be exactly one response variable")
data.env <- environment(Fo)
if(!is.null(data)){
if(is.list(data)){ # list (including data.frame)
data.env <- list2env(data, parent = data.env)
environment(Fo) <- data.env
# put the response in y.
# TODO: I had forgotten that I do this! Is it a good idea?
data <- model.part(Fo, data = data, lhs = 1)
data <- as.matrix(data) # TODO: should this be as.data.frame?
}else if(is.environment(data))
stop("'data' from class 'environment' not supported")
else{ # data is the time series # TODO: consider multivariate
data.env <- new.env(parent = data.env)
environment(Fo) <- data.env
data.env[[yname]] <- data
mf0 <- model.frame(Fo, lhs = 1, rhs = FALSE)
data <- model.part(Fo, data = mf0, lhs = 1)
data <- as.matrix(data) # TODO: should this be as.data.frame?
}
}else{
# !!! TODO: do this for the above as well !!!
#
# 15/01/2018 - na.action = NULL, otherwise NA's are dropped
# mf0 <- model.frame(Fo, lhs = 1, rhs = FALSE)
mf0 <- model.frame(Fo, lhs = 1, rhs = FALSE, na.action = NULL)
data <- model.part(Fo, data = mf0, lhs = 1)
data <- as.matrix(data) # TODO: should this be as.data.frame?
}
stopifnot(identical(environment(Fo), data.env))
# get nobs - there should be a better way, since
# calling model frame requires that the
# correct response variable is visible. I
# work only with the lhs, so the other
# variables are not needed at this stage
# (as it have to - cs, ets., are set
# using nobs further below).
Fo.0 <- formula(Fo, lhs = 1, rhs = FALSE)
## 15/01/2018 - bug fix: this drops NA's!
## nobs <- nrow(model.frame(Fo.0)) # TODO: model.frame() needs the environment to be set up!
# any other wy to find nobs?
# DONE: Moving stuff seting the environment before this line
# but leaving this note here in case I forget.
nobs <- nrow(model.frame(Fo.0, na.action = NULL))
#browser()
## NOTE: Do not reuse Fo.0 - its environment gets out of sync with Fo.0o after the following.
## (but the data are already there, so maybe no problem.
# insert an environment for the trend stuff
data.env <- new.env(hash = FALSE, parent = environment(Fo))
data.env$t <- seq_len(nobs)
Fo.trend <- formula(Fo, rhs = 1)
Fo.trend.rhs <- formula(Fo, lhs = FALSE, rhs = 1)
trmake <- NULL
if(!is.empty.model(Fo.trend))
# lhs is needed if there is only intercept, otherwise nobs is unknown
# trmake <- trendMaker(Fo.trend.rhs, time = 1:nobs)
trmake <- trendMaker(Fo.trend, time = 1:nobs)
Fo.sarima <- formula(Fo, rhs = 2) # TODO: check for empty
udelta <- delta <- theta <- phi <- list() # initialise variables, will be updated by toPolyCoef()
molist <- SARIMA(Fo.sarima)
lapply(molist, toPolyCoef)
## renaming Xreg to regx
regxmake <- NULL
if(length(Fo)[2] > 2){
Fo.regx <- formula(Fo, rhs = 3)
if(!is.empty.model(Fo.regx))
regxmake <- trendMaker(Fo.regx, time = 1:nobs)
}else{
Fo.regx <- NULL # for clarity; currently not used
}
#browser()
res <- sarimat(data, phi, theta, delta, udelta, trmake = trmake, regxmake = regxmake,
lik.method = lik.method, use.symm = use.symmetry, SSinit = SSinit)
res$call <- match.call() # replace the call to sarimat() with the call to sarima()
res$formula <- model
res$internal$Fo.xreg <- Fo.trend
res$internal$Fo.sarima <- formula(Fo, lhs = FALSE, rhs = 2) # Fo.sarima
res$internal$Fo.regx <- Fo.regx
#browser()
res
}
sarimat <- function(y, phi, theta, delta, udelta, trmake = NULL, regxmake = NULL, lik.method,
use.symm = FALSE, SSinit = c("Rossignol2011", "gnb", "Gardner1980")){
## 2018-10-12
## SSinit_Q0 <- "Rossignol2011"
## SSinit_Q0 <- "Gardner1980"
SSinit_Q0 <- match.arg(SSinit)
makeArimaInternal <- if(SSinit_Q0 == "gnb")
makeArimaGnb
else
makeARIMA
###### JAMIE: changed 23/08/2018 to include a symmetric approach to U(z)
update_params <- function(par, use.symm = FALSE){
## flat_par contains the parameters as seen by optim,
## except that optim only gets the nonfixed parameters
flat_par[nonfixed] <<- par
## the calculation below depend on the "real" parameters,
## so we tanh transform if necessary to get parcor's, then transform the parcor's
## to coefficients.
## this assumes that the sarima parameters are before any others
## (which is the case, but this is are minder not to change it!)
## using integer indices rather than logical ones to avoid extending
## flags.atanh_sarima to the full length of flat_par.
##
## TODO: In some cases this transformation may make sense for other parameters too.
## But let's sort the sarima parameters first
if(tanh.flag){ ## tanh() parameters to get parcor's
wrk <- flat_par
## fixed parameters declared to need transform, will be transformed to and from.
## TODO: fix this
wrk[ind.atanh_sarima] <- tanh(flat_par[ind.atanh_sarima])
par <- wrk[nonfixed]
}
newpars <- relist(par, par.skeleton)
mapply(function(x, newpar){
x$parcor[x$dyn.index] <- newpar
x$a[-1] <- - pacf2Ar(x$parcor[-1]) # 2022-02-14 was: FitAR::PacfToAR(x$parcor[-1])
}, e_phi, newpars[["phi"]], SIMPLIFY = FALSE)
mapply(function(x, newpar){
x$parcor[x$dyn.index] <- newpar
x$a[-1] <- pacf2Ar(x$parcor[-1]) # 2022-02-14 was: FitAR::PacfToAR(x$parcor[-1])
}, e_theta, newpars[["theta"]], SIMPLIFY = FALSE)
if(use.symm){
mapply(function(x, newpar){
x$parcor_symm[x$dyn.index_symm] <- newpar
x$parcor[x$polyInd] <- x$parcor_symm
x$a[-1] <- - pacf2Ar(x$parcor[-1]) # 2022-02-14 was: FitAR::PacfToAR(x$parcor[-1])
x$coef_symm <- x$a[x$paramInd]
}, e_udelta, newpars[["udelta"]], SIMPLIFY = FALSE)
}else{
mapply(function(x, newpar){
x$parcor[x$dyn.index] <- newpar
x$a[-1] <- - pacf2Ar(x$parcor[-1]) # 2022-02-14 was: FitAR::PacfToAR(x$parcor[-1])
}, e_udelta, newpars[["udelta"]], SIMPLIFY = FALSE)
}
Par_treg <<- if(treg.flag) unlist(newpars[["xreg"]])
else numeric(0)
Par_const <<- if(X.flag) unlist(newpars[["regx"]])
else numeric(0)
NULL
}
makePhiTheta <- function(){
Phi <<- - coef(sprod(phi_poly))[-1]
if(length(Phi) != Phi_length)
Phi <<- c(Phi, numeric(Phi_length - length(Phi)))
Theta <<- coef(sprod(theta_poly))[-1]
if(length(Theta) != Theta_length)
Theta <<- c(Theta, numeric(Theta_length - length(Theta)))
Phi <<- - coef(sprod(phi_poly))[-1]
Delta_prod <<- delta_prod * sprod(udelta_poly)
Delta <<- - coef(Delta_prod)[-1]
if(length(Delta) != Delta_length)
Delta <<- c(Delta, numeric(Delta_length - length(Delta)))
}
makeArimaConst <- function(const, ...){
## see Durbin/Koopman, p. 61
## only const for now, since otherwise Z would be time dependent.
mod <- makeArimaInternal(...)
mod$Z <- c(1, mod$Z)
mod$a <- c(const, mod$a)
## not Matrix::bdiag() for a reason
## - package currently doesn't import Matrix
## - needs to be able to handle e.g. a 0-row matrix
mod$P <- rbind(c(const, numeric(ncol(mod$P))), cbind(0, mod$P))
mod$Pn <- rbind(c(0, numeric(ncol(mod$Pn))), cbind(0, mod$Pn))
mod$T <- rbind(c(1, numeric(ncol(mod$T))), cbind(0, mod$T))
## V = RQR', diag_bind(0, V), i.e. non-random intercept
mod$V <- rbind(c(0, numeric(ncol(mod$V))),
cbind(0, mod$V))
mod
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
loglik <- function(par, use.symm = FALSE){
update_params(par, use.symm)
makePhiTheta() # expand polynomials
if(X.flag){
model_from_makeARIMA <<- makeArimaConst(Par_const, Phi, Theta, Delta, SSinit = SSinit_Q0)
}else
model_from_makeARIMA <<- makeArimaInternal(Phi, Theta, Delta, SSinit = SSinit_Q0)
if(treg.flag)
yc <<- y - tregmat %*% Par_treg
#browser()
fromKalman <<- KalmanLike(yc, model_from_makeARIMA)
fromKalman$Lik
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
loglik_sarima <- function(par, use.symm = FALSE){
update_params(par, use.symm)
makePhiTheta() # expand polynomials
if(X.flag){
model_from_makeARIMA <<- makeArimaConst(Par_const, Phi, Theta, Delta, SSinit = SSinit_Q0)
}else
model_from_makeARIMA <<- makeArimaInternal(Phi, Theta, Delta, SSinit = SSinit_Q0)
if(treg.flag)
yc <<- y - tregmat %*% Par_treg
# fromKalman <<- KalmanLike(yc, model_from_makeARIMA)
fromKalman <<- sarima_KalmanLike(yc, model_from_makeARIMA)
#cat("Phi = ", Phi, " Theta = ", Theta, "\n")
#browser()
fromKalman$Lik
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
ss_sarima <- function(par, use.symm = FALSE){
update_params(par, use.symm)
makePhiTheta() # expand polynomials
if(X.flag){
model_from_makeARIMA <<- makeArimaConst(Par_const, Phi, Theta, Delta, SSinit = SSinit_Q0)
}else
model_from_makeARIMA <<- makeArimaInternal(Phi, Theta, Delta, SSinit = SSinit_Q0)
if(treg.flag)
yc <<- y - tregmat %*% Par_treg
# fromKalman <<- KalmanLike(yc, model_from_makeARIMA)
fromKalman <<- sarima_KalmanLike(yc, model_from_makeARIMA)
fromKalman$css
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
loglik_fkf <- function(par, use.symm = FALSE){
npar.fkf <- length(par)
## @jamie TODO: You set it to this:
## update_params(par, use.symm)
## However, I leave the original, since it doesn't pass all parameters:
## TODO: check if this is oversight.
##
## @georgi (10/10/18) This is as it should be.
##
update_params(par[- npar.fkf])
sigma <<- exp(par[npar.fkf])
makePhiTheta() # expand polynomials
## TODO: check X.flag, etc.
ct <- t(tregmat %*% Par_treg)
dt <- t(regxmat %*% Par_const)
fkf.model <<- xarimaxss(ar = Phi, ma = Theta, sigma = sigma,
xreg = ct, regx = dt, delta = Delta)
fromKalman <<- FKF::fkf(a0 = fkf.model$a0, P0 = fkf.model$P0,
dt = fkf.model$dt, ct = fkf.model$ct,
Tt = fkf.model$Tt, Zt = fkf.model$Zt,
HHt = fkf.model$HHt, GGt = fkf.model$GGt,
yt = matrix(y, nrow = 1))
- fromKalman$logLik # TODO: - $logLik, while in "base" it is $Lik !!!
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
loglik_kfas <- function(par, use.symm = FALSE){
npar.fkf <- length(par) # last param is log(sigma2) (In FKF it is log(sigma))
## @jamie TODO: leaving as is,see note above
## @georgi (10/10/18): This is as it should be - as above.
update_params(par[- npar.fkf])
# ## TODO: correct the hessian later
# beta <- tanh(par[- npar.fkf])
# cat("beta = ", beta, "\n")
#
# update_params_kfas(beta)
## TODO: is it better to use sigma as parameter for optim()?
## (better scaling?)
sigma2 <<- exp(par[npar.fkf]) # TODO: check - is this sigma2 or sigma squared?
makePhiTheta() # expand polynomials
# if(X.flag){
# model_from_makeARIMA <<- makeArimaConst(Par_const, Phi, Theta, Delta)
# }else
# model_from_makeARIMA <<- makeARIMA(Phi, Theta, Delta)
# if(treg.flag)
# yc <<- y - tregmat %*% Par_treg
#
# fromKalman <<- KalmanLike(yc, model_from_makeARIMA)
#
# fromKalman$Lik
## TODO: integrated, trends etc,
# ct <- t(tregmat %*% Par_treg)
# dt <- t(regxmat %*% Par_const)
# fkf.model <<- xarimaxss(ar = Phi, ma = Theta, sigma = sigma,
# xreg = ct, regx = dt, delta = Delta)
# from makeUpdateFun_arma0() and KFAS examples
# TODO: only ARMA w.o. intercept for now.
wrk <- try(KFAS::SSMarima(Phi, Theta, Q = sigma2), silent = TRUE)
if(inherits(wrk, "try-error")){
# cat("1e100 !!\n")
return(1e100)
}
kfas_model["T", "arima"] <<- wrk$T
kfas_model["R", "arima"] <<- wrk$R
kfas_model["P1", "arima"] <<- wrk$P1
kfas_model["Q", "arima"] <<- wrk$Q
fromKalman <<- list(logLik = logLik(kfas_model))
- fromKalman$logLik # TODO: - $logLik, while in "base" it is $Lik !!!
}
poly_param <- function(x, fixed = attr(x, "fixed"), nseasons = attr(x, "nseasons"),
d = attr(x, "d", exact = TRUE),
atanh.tr = attr(x, "atanh.tr", exact = TRUE) ){
coef <- if(!is.null(attr(x, "sign")) && attr(x, "sign") == "-")
c(1, - x)
else
c(1, x)
# polynom() drops trailing zeroes and doesn't accept NA's.
# to protect against these, create a polynomial of 1's first
res <- polynom(rep(1, length(coef)))
## TODO: krapka - replace NA's with zeroes
if(anyNA(coef))
coef[is.na(coef)] <- 0
e <- environment(res)
e$a <- coef
e$fixed <- if(is.null(fixed))
c(TRUE, rep(FALSE, length(coef) - 1))
else
c(TRUE, fixed)
e$dyn.index <- which(!e$fixed)
## 18/01/2017 - support for atanh transform
## TODO: maybe different defaults for AR and MA? but here it is not known if
## the parameter is AR or MA.
## It is best for the default here to be FALSE since otherwise the caller
## will have to take care of this even for model elements for which this
## doesn't make sense. In any case, the default here is not that important
## since sarima() may set these.
e$atanh.tr <- if(is.null(atanh.tr))
c(FALSE, rep(FALSE, length(coef) - 1))
else
c(FALSE, atanh.tr)
## TODO: for now there seems no need for analog of dyn.index (see above)
if(!is.null(nseasons)){
## TODO: check for whole number
if(is.numeric(nseasons) && length(nseasons) == 1 && nseasons > 0){
e$nseasons <- nseasons
}else
stop("'nseasons' must be a positive integer or NULL")
}
if(!is.null(d) && d > 1){
e$d <- d
}
if(!is.null(dnam <- attr(x, "dispname")))
e$dispname <- dnam
res
}
sprod <- function(pp){ # pp is short for parameter polynomial
f <- function(poly){
s <- environment(poly)$nseasons
if(is.null(s))
poly
else
poly(poly_x^s)
}
prod(as.polylist(lapply(pp, f)))
}
u_sprod <- function(pp){ # pp is short for parameter polynomial
f <- function(poly){
s <- environment(poly)$nseasons
if(!is.null(s))
poly <- poly(poly_x^s)
d <- environment(poly)$d
if(is.null(d)) poly else poly^d
}
prod(as.polylist(lapply(pp, f)))
}
sdegree <- function(pp){
e <- environment(pp)
s <- e$nseasons
degree <- length(e$a) - 1
if(!is.null(s))
degree <- degree * s
if(!is.null(e$d))
degree <- degree * e$d
degree
}
## TODO: this is krapka for naming multple ar() terms, etc.
poly_param_plus <- function(parlist){
no <- 0
sno <- 0
res <- vector(length(parlist), mode = "list")
for(i in seq_along(parlist)){
wrk <- poly_param(parlist[[i]])
e <- environment(wrk)
if(is.null(e$nseasons)){
e$name.counter <- no
no <- no + 1
}else{
e$name.counter <- sno
sno <- sno + 1
}
res[[i]] <- wrk
}
as.polylist(res)
}
# phi, theta and delta are used only here.
# TODO: maybe the they can come as polynomials to sarimat() ?
phi_poly <- poly_param_plus(phi)
theta_poly <- poly_param_plus(theta)
delta_poly <- poly_param_plus(delta)
udelta_poly <- poly_param_plus(udelta)
# Delta is fixed, so needs to be computed only once
## TODO: (1-x)^d is kept as 1-x and a 'd' is a sepaarate varible.
## the element of delta_poly corresponding to this will print as 1 - x.
## u_sprod(), print.Sarima(), summary.Sarima() know this but ing eneral this is very fragile.
delta_prod <- u_sprod(delta_poly)
udelta_prod <- sprod(udelta_poly) # NOTE: this is not u_sprod(udelta_poly)
Par_const <- Par_treg <- Phi <- Theta <- numeric(0)
yc <- y
sigma <- sd(y) # for fkf
sigma2 <- var(y) # for kfas
kfas_model <- fkf.model <- fromKalman <- model_from_makeARIMA <- NULL
Delta_prod <- delta_prod * udelta_prod
Delta <- - coef(Delta_prod)[-1]
e_phi <- lapply(phi_poly, environment)
e_theta <- lapply(theta_poly, environment)
e_delta <- lapply(delta_poly, environment)
e_udelta <- lapply(udelta_poly, environment)
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - add extra arguments to the environment to help this approach
if(use.symm){
.nonstationaryEdit <- function(e){
## Parcors will be symmetric inside the extreme parcors
## so we only want half the parameters
len <- length(e$a)
midpoint <- ceiling(len/2 - 1)
mid_seq <- seq_len(midpoint)
## fixed and atanh.tr, accounting for symmetry
coef_a <- e$a[-c(1,len)]
fixed_a <- e$fixed[-c(1,len)]
atanh_a <- e$atanh.tr[-c(1,len)]
fixed_symm <- fixed_a[mid_seq] | rev(fixed_a)[mid_seq]
atanh_symm <- atanh_a[mid_seq] | rev(atanh_a)[mid_seq]
## where to find b parameters
paramInd <- mid_seq + 1
## starting coefficients
coef_symm <- e$a[-c(1,len)][mid_seq]
## if any value is fixed, make sure we take it appropriately
if(any(fixed_symm)){
fix_a <- which(fixed_a)
for(i in seq_along(fix_a)){
if(fix_a[i] > midpoint){
fix_symm <- len - 2 - fix_a[i]
paramInd[fix_symm] <- fix_a[i] + 1
}else{
fix_symm <- fix_a[i]
}
coef_symm[fix_symm] <- coef_a[fix_a[i]]
}
}
## add to environment
e$fixed_symm <- fixed_symm
e$atanh.tr_symm <- atanh_symm
e$dyn.index_symm <- which(!fixed_symm)
e$paramInd <- paramInd
e$coef_symm <- coef_symm
## where b parameters can be found
e$polyInd <- c(seq_along(coef_symm) + 1,
len - seq_along(coef_symm))
NULL
}
lapply(e_udelta, .nonstationaryEdit)
}
e_phi_theta <- c(e_phi, e_theta)
e_phi_theta_delta <- c(e_phi_theta, e_delta)
e_phi_theta_udelta <- c(e_phi_theta, e_udelta)
Phi_length <- if(length(phi_poly) > 0) sum(sapply(phi_poly, sdegree)) else 0
Theta_length <- if(length(theta_poly) > 0) sum(sapply(theta_poly, sdegree)) else 0
Udelta_length <- if(length(udelta_poly) > 0) sum(sapply(udelta_poly, sdegree)) else 0
delta_length <- if(length(delta_poly) > 0) sum(sapply(delta_poly, sdegree)) else 0
Delta_length <- delta_length + Udelta_length
prepare_parcor <- function(envir){
if(is.null(envir$parcor)){
envir$parcor <- c(1, ar2Pacf(- envir$a[-1]))
}else if(length(envir$parcor) != length(envir$a))
stop("lengths of coef and parcor are not equal")
if(anyNA(envir$parcor)) # todo: give a warning in this case?
envir$parcor[is.na(envir$parcor)] <- 0
NULL
}
# set initial values for parcor's (transformed params)
# also: "parcor's don't make sense for delta!
lapply(e_phi_theta, prepare_parcor)
lapply(e_udelta, function(e) e$parcor <- e$a) # TODO: $parcor should be set outright, this is a hack
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - update the parameters that will be driving the equation
if(use.symm) lapply(e_udelta, function(e) e$parcor_symm <- e$parcor[e$paramInd])
# TODO: this is just wrong, but not used:
lapply(e_delta, function(e) e$parcor <- e$a) # TODO:delta !!!!!
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - The symmetric approach uses fewer parameters, so the length
###### of variables in the functions defined at the start needs to
###### be the correct size, etc.
if(use.symm){
flags.atanh_phi_theta <- lapply(e_phi_theta, function(x)
if(isTRUE(x$operator)) logical(0) else x$atanh.tr)
flags.atanh_udelta <- lapply(e_udelta, function(x)
if(isTRUE(x$operator)) logical(0) else x$atanh.tr_symm)
flags.atanh_sarima <- unlist(c(flags.atanh_phi_theta, flags.atanh_udelta))
if(length(flags.atanh_sarima) == 0) ## in case it is list()
flags.atanh_sarima <- logical(0)
ind.atanh_sarima <- which(flags.atanh_sarima)
tanh.flag <- any(flags.atanh_sarima)
Par_phi_theta <- lapply(e_phi_theta, function(x)
if(isTRUE(x$operator)) numeric(0) else x$parcor)
Par_udelta <- lapply(e_udelta, function(x)
if(isTRUE(x$operator)) numeric(0) else x$parcor_symm)
Par_sarima <- unlist(c(Par_phi_theta, Par_udelta))
fixed_phi_theta <- lapply(e_phi_theta, function(x)
if(isTRUE(x$operator)) logical(0) else x$fixed)
fixed_udelta <- lapply(e_udelta, function(x)
if(isTRUE(x$operator)) logical(0) else x$fixed_symm)
fixed <- unlist(c(fixed_phi_theta, fixed_udelta))
}else{
flags.atanh_sarima <- unlist(lapply(e_phi_theta_udelta,
function(x)
if(isTRUE(x$operator)) logical(0) else x$atanh.tr
))
if(length(flags.atanh_sarima) == 0) ## in case it is list()
flags.atanh_sarima <- logical(0)
ind.atanh_sarima <- which(flags.atanh_sarima)
tanh.flag <- any(flags.atanh_sarima)
Par_sarima <- unlist(lapply(e_phi_theta_udelta,
function(x)
if(isTRUE(x$operator)) numeric(0) else x$parcor
))
## @jamie TODO: !!!
## You have changed the logic even n the non-sym case, here and at other places.
## It is extremely dangerous to do multiple unrelated sweeping changes
## at the same time, especcially without proper testing.
## e.g. this assignement has been mmoved from further below here, silently.
fixed <- unlist(lapply(e_phi_theta_udelta,
function(x)
if(isTRUE(x$operator)) logical(0) else x$fixed
))
}
stopifnot(length(Par_sarima) == length(flags.atanh_sarima))
if(tanh.flag)
Par_sarima[flags.atanh_sarima] <- atanh(Par_sarima[flags.atanh_sarima])
X.flag <- !is.null(regxmake)
treg.flag <- !is.null(trmake)
## now prepare the treg parameters
## time regression
## set up xreg
if(treg.flag){
tregmat <- trmake()
if(length(Delta) > 0){
treg.temp <- trmake((1 - length(Delta)) : nrow(tregmat))
treg.temp <- cbind(c(rep(NA_real_, length(Delta)), y), treg.temp)
# TODO: use filter()
## start from the last row to do it in-place
for(i in nrow(treg.temp):(length(Delta) + 1)){
for(j in 1:length(Delta)){
treg.temp[i, ] <- treg.temp[i, ] - Delta[j] * treg.temp[i - j, ]
}
}
## drop the added rows
treg.temp <- treg.temp[-(1:length(Delta)), ]
ydiff <- treg.temp[-(1:length(Delta)), 1]
tregdiff <- treg.temp[-(1:length(Delta)), -1]
fitinit <- lm(ydiff ~ tregdiff - 1, na.action = na.omit)
Par_treg <- coef(fitinit)
}else{
# Tova nyama da e vyarno ako ima regx:
# if(ncol(tregmat) == 1 && all(tregmat == 1)){
# colnames(tregmat) <- "mean"
# }
fitinit <- lm(y ~ tregmat - 1, na.action = na.omit)
Par_treg <- coef(fitinit)
}
}else{#treg.flag = FALSE
tregmat <- matrix(0, nrow = length(y), ncol = 0)
Par_treg <- numeric(0)
}
## set up regx
if(X.flag){
regxmat <- regxmake()
fitinit <- lm(y ~ regxmat - 1, na.action = na.omit)
Par_const <- coef(fitinit)
}else{
regxmat <- matrix(0, nrow = length(y), ncol = 0)
Par_const <- numeric(0)
}
# note: flat_par includes the coeffcients of zero power of the polynomials
flat_par <- c(Par_sarima, Par_treg, Par_const)
# TODO: for now no fixed treg parameters
# use 'offset' if possible?
Par_treg_fixed <- rep(FALSE, length(Par_treg))
Par_const_fixed <- rep(FALSE, length(Par_const))
fixed <- c(fixed, Par_treg_fixed, Par_const_fixed)
## note: flat_par includes the coeffcients of zero power of the polynomials;
## so, flat_par[nonfixed] includes them, too!
nonfixed <- !fixed
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - par.skeleton[["udelta]] will have fewer dynamic variables
if(use.symm){
par.skeleton <- list(phi = lapply(e_phi, function(e) e$dyn.index),
theta = lapply(e_theta, function(e) e$dyn.index),
udelta = lapply(e_udelta, function(e) e$dyn.index_symm),
# delta = lapply(e_delta, function(e) e$dyn.index),
xreg = list(seq_along(Par_treg)),
regx = list(seq_along(Par_const)))
}else{
par.skeleton <- list(phi = lapply(e_phi, function(e) e$dyn.index),
theta = lapply(e_theta, function(e) e$dyn.index),
udelta = lapply(e_udelta, function(e) e$dyn.index),
# delta = lapply(e_delta, function(e) e$dyn.index),
xreg = list(seq_along(Par_treg)),
regx = list(seq_along(Par_const)))
}
## optimise
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - each method now includes an argument to use a symmetric approach
switch(lik.method,
base = {
loglik(flat_par[nonfixed], use.symm = use.symm)
res_opt <- optim(flat_par[nonfixed], loglik, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
# update the quantities with the result from the optim:
# (the last call of loglik() by optim() may not be the final result)
loglik(res_opt$par, use.symm = use.symm) # identical(res_opt$value, fromKalman$Lik) should be TRUE
},
sarima = {
loglik_sarima(flat_par[nonfixed], use.symm = use.symm)
res_opt <- optim(flat_par[nonfixed], loglik_sarima, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
if(tanh.flag){
## lazy: call optim with the untransformed parameters with zero maxit,
## simply to get the hessian for the untransformed parameters.
## todo: however the current setup makes this difficult
## save some vars
tanh.flag <- FALSE
sav_flat_par <- flat_par
## tanh() to get parcor's
flat_par[ind.atanh_sarima] <- tanh(flat_par[ind.atanh_sarima])
## this is only needed for the hessian
res_opt_alt <- optim(flat_par[nonfixed], loglik_sarima, use.symm = use.symm,
method = "BFGS", hessian = TRUE,
control = list(maxit = 0L))
flat_par <- sav_flat_par # restore saved vars
tanh.flag <- TRUE
}
# update the quantities with the result from the optim:
# (the last call of loglik() by optim() may not be the final result)
loglik_sarima(res_opt$par, use.symm = use.symm) # identical(res_opt$value, fromKalman$Lik) should be TRUE
},
css = {
ss_sarima(flat_par[nonfixed], use.symm = use.symm)
res_opt <- optim(flat_par[nonfixed], ss_sarima, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
if(tanh.flag){
## save some vars
tanh.flag <- FALSE
sav_flat_par <- flat_par
## tanh() to get parcor's
flat_par[ind.atanh_sarima] <- tanh(flat_par[ind.atanh_sarima])
## this is only needed for the hessian
res_opt_alt <- optim(flat_par[nonfixed], ss_sarima, use.symm = use.symm,
method = "BFGS", hessian = TRUE,
control = list(maxit = 0L))
flat_par <- sav_flat_par # restore saved vars
tanh.flag <- TRUE
}
ss_sarima(res_opt$par, use.symm = use.symm)
## 10/10/18 method changed to 'sarima' to avoid duplicating code later
lik.method <- "sarima"
},
fkf = {
loglik_fkf(c(flat_par[nonfixed], log(sigma)), use.symm = use.symm)
res_opt <- optim(c(flat_par[nonfixed], log(sigma)),
loglik_fkf, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
loglik_fkf(res_opt$par, use.symm = use.symm)
},
kfas = {
# KFAS::SSMarima() returns a list
# kfas_model <- try(KFAS::SSMarima(Phi, Theta, Q = sigma2), silent = TRUE)
## loglik_kfas() doesn't recreate the full model,
## so create it in advance (could add if(is.null(kfas_model)) ... to loglik_kfas)
update_params(flat_par[nonfixed], use.symm = use.symm)
makePhiTheta() # expand polynomials
## 2022-01-19 moved KFAS to Suggests.
##
## however, then SSMarima is not visible and in a model formula can't
## be qualified with 'KFAS::'. So, copying it here.
## TODO: move this assignment elsewhere (here it is called every
## time the likelihood is computed, the cost is negligible though)
SSMarima <- KFAS::SSMarima
kfas_model <- KFAS::SSModel(y ~ -1 + SSMarima(Phi, Theta, Q = sigma2), H = 0)
# kfas_init <- c(atanh(flat_par[nonfixed]), log(sigma2))
kfas_init <- c(flat_par[nonfixed], log(sigma2))
loglik_kfas(kfas_init, use.symm = use.symm)
res_opt <- optim(kfas_init, loglik_kfas, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
loglik_kfas(res_opt$par, use.symm = use.symm)
},
arima = {
## this methods calls arima()
if(length(Par_const) > 0)
stop("method 'arima' can handle 'xreg', but not 'regx'")
## TODO: finish this
stop("method arima is unfinished")
},
stop("invalid 'lik.method'")
)
fcoef <- function(x, minus = FALSE){
res <- x$a[-1]
if(minus)
res <- - res
if(!is.null(nam <- x$dispname)){
suffix <- letters[x$name.counter] # should be ok for NULL and 0
names(res) <- paste0(nam, seq_along(res), suffix)
}
res
}
coef_phi <- lapply(e_phi, function(x) fcoef(x, TRUE))
fixed_phi <- lapply(e_phi, function(x) x$fixed[-1])
coef_theta <- lapply(e_theta, fcoef)
fixed_theta <- lapply(e_theta, function(x) x$fixed[-1])
coef_delta <- lapply(e_delta, function(x) fcoef(x, TRUE))
fixed_delta <- lapply(e_delta, function(x) x$fixed[-1])
coef_udelta <- lapply(e_udelta, function(x) fcoef(x, TRUE))
fixed_udelta <- lapply(e_udelta, function(x) x$fixed[-1])
all.coef <- unlist(c(coef_phi, coef_theta, coef_udelta, coef_delta, Par_treg, Par_const))
all.fixed <- unlist(c(fixed_phi, fixed_theta, fixed_udelta, fixed_delta, Par_treg_fixed, Par_const_fixed))
coef.nonfixed <- unlist(all.coef)[!all.fixed]
# TODO: think about better names, can clean them up here
if(length(Par_treg) > 0)
names(Par_treg) <- colnames(tregmat)
if(length(Par_const) > 0)
names(Par_const) <- colnames(regxmat)
# no delta here, but all this needs further thought
# TODO: need to add numbers to c("ar", "ar"), etc
coef <- unlist(c(coef_phi, coef_theta, coef_udelta, Par_treg, Par_const))
fixed_coef <- unlist(c(fixed_phi, fixed_theta, fixed_udelta, Par_treg_fixed, Par_const_fixed))
## TODO: Tova e krapka, ako e list() ili NULL , !fixed_coef po-dolu dava greshka.
if(length(fixed_coef) == 0)
fixed_coef <- logical(0)
n.coef <- c(arma = length(unlist(c(coef_phi, coef_theta))),
# TODO: reinstate but change predict.Sarima to use indexing with names
uar = length(unlist(coef_udelta)),
xreg = length(Par_treg),
regx = length(Par_const))
mask <- !fixed_coef
## n.cond in arima is only for method css, maybe should drop it here.
ncond <- 0 # max(length(Phi), length(Theta) + 1) # TODO: Fix this!
# in arima() n.used is:
# n.used <- sum(!is.na(x)) - length(Delta)
# (with additional correction for missing values in xreg).
# TODO: do the correction for xreg.
notna <- !is.na(y)
n.used <- if(length(Delta) == 0)
sum(notna)
else if(all(!is.na(y[1:length(Delta)])))
sum(notna) - length(Delta)
else{ # NA's among the first length(Delta) values - additional loss
# if there is a stretch of length(Delta) non-NA's all subsequent obs. are ok
# but the exact number of non-ok depends also on zero coefficients in Delta
# (e.g. for 1 - B^12)
## TODO: remove this or further modify at least for ss.method = "sarima"
available <- !is.na(y)
statable <- rep(FALSE, length(y))
ddind <- which(Delta != 0)
okseq <- 0
for(t in (length(Delta) + 1):length(y)){
statable[t] <- all(statable[t - ddind] | available[t - ddind])
if(statable[t])
okseq <- okseq + 1
else
okseq <- 0
if(okseq >= length(Delta)){
if(t < length(y))
statable[(t + 1):length(y)] <- TRUE
break
}
}
## !!! TODO: statable may be useful to choose which residuals to use etc.
sum(statable)
}
switch(lik.method,
base = {
# arima(): value <- 2 * n.used * res$value + n.used + n.used * log(2 * pi)
# aic <- if(method != "css") value + 2*sum(mask) + 2 else NA
# KalmanLike() returns minus log-likelihood, upto a scaling and addition of constants
# (not fully documented, the following is a guess from trial and error and regression)
# equivalently, can replace fromKalman$Lik with res_opt$value
value_loglik <- - n.used * fromKalman$Lik - 0.5 * n.used * log(2 * pi) - 0.5 * n.used
# I don't fully understand argument 'nit', but setting it to zero seems
# to be right for prediction (setting it to -1 gives different results), compare:
# fromKalmanRun <- KalmanRun(y, model_from_makeARIMA, -1, update = TRUE)
# fromKalmanRun <- KalmanRun(y, model_from_makeARIMA, 0, update = TRUE)
fromKalmanRun <- KalmanRun(yc, model_from_makeARIMA, 0, update = TRUE)
resid <- fromKalmanRun$resid
### from jamie: these are standardised residuals!
# TODO: sigma2 needs more care.
# This reproduces (it seems) arima()'s sigma2 in the absence of NA's
ndrop_sigma2 <- length(Delta)
sigma2 <- if(ndrop_sigma2 == 0)
mean(resid^2, na.rm = TRUE)
else
mean(resid[-seq_len(ndrop_sigma2)]^2, na.rm = TRUE)
sigma <- sqrt(sigma2)
## should be ready for prediction (minus exogeneous vars)
ss.model <- attr(fromKalmanRun, "mod") # model_from_makeARIMA
coef_hessian <- res_opt$hessian
#browser()
},
sarima = {
## for now same as "base" but calling my likelihood code
value_loglik <- - n.used * fromKalman$Lik - 0.5 * n.used * log(2 * pi) - 0.5 * n.used
# I don't fully understand argument 'nit', but setting it to zero seems
# to be right for prediction (setting it to -1 gives different results), compare:
# fromKalmanRun <- KalmanRun(y, model_from_makeARIMA, -1, update = TRUE)
# fromKalmanRun <- KalmanRun(y, model_from_makeARIMA, 0, update = TRUE)
## TODO: this needs changing, but lets first sort out the estimation
## The C++ function uniKalmanLikelihood0b needs amending to return states and
## the updated model
## fromKalmanRun <- sarima_KalmanRun(yc, model_from_makeARIMA, 0, update = TRUE)
## (2018-08-23) TODO:
## @jamie: In such cases please put a note has happened here.
## Is this a follow up to the aboe note, if so to what extent, etc.
## I assume that you have tested the equivalence of these two in the context.
## 2018-08-23 was:
## fromKalmanRun <- KalmanRun(yc, model_from_makeARIMA, 0, update = TRUE)
fromKalmanRun <- sarima_KalmanRun(yc, model_from_makeARIMA, 0, update = TRUE)
resid <- fromKalmanRun$resid
# TODO: sigma2 needs more care.
# This reproduces (it seems) arima()'s sigma2 in the absence of NA's
ndrop_sigma2 <- length(Delta)
sigma2 <- if(ndrop_sigma2 == 0)
mean(resid^2, na.rm = TRUE)
else
mean(resid[-seq_len(ndrop_sigma2)]^2, na.rm = TRUE)
sigma <- sqrt(sigma2)
## should be ready for prediction (minus exogeneous vars)
ss.model <- attr(fromKalmanRun, "mod") # model_from_makeARIMA
## TODO: it may be better not to call again optim
## but to transform the hessian as in fitArma0Model()
coef_hessian <- if(tanh.flag)
res_opt_alt$hessian
else
res_opt$hessian
},
fkf = {
value_loglik <- fromKalman$logLik
## fkf puts other useful stuff in the answers
resid <- as.vector(fromKalman$vt)
sigma2 <- sigma^2
ss.model <- fkf.model
nrhess <- dim(res_opt$hessian)[1] # hess is square matrix
coef_hessian <- res_opt$hessian[- nrhess, - nrhess] # drop logsigma
},
kfas = {
value_loglik <- fromKalman$logLik
## fkf puts other useful stuff in the answers
resid <- as.vector(fromKalman$vt) # TODO: this is incorrect!
ss.model <- kfas_model
nrhess <- dim(res_opt$hessian)[1] # hess is square matrix
coef_hessian <- res_opt$hessian[- nrhess, - nrhess] # drop logsigma
},
arima = {
stop("Should not get this far, 'arima' method not supported yet.")
},
stop("'lik.method' not found.")
)
aic <- - 2 * value_loglik + 2 * sum(nonfixed) + 2
transformed_par <- list(lapply(e_phi, function(x) x$parcor),
lapply(e_theta, function(x) x$parcor),
lapply(e_udelta, function(x) x$parcor),
lapply(e_delta, function(x) x$parcor) )
phipar_only <- unlist(lapply(e_phi, function(x) x$a[-1]))
phipar_fixed <- unlist(lapply(e_phi, function(x) x$fixed[-1]))
thetapar_only <- unlist(lapply(e_theta, function(x) x$a[-1]))
thetapar_fixed <- unlist(lapply(e_theta, function(x) x$fixed[-1]))
# Jphi, Jtheta are lists of matrices
# transformed par are the parcor's
Jphi <- lapply(transformed_par[[1]], function(x) pacf2ArWithJacobian(x[-1])$J )
# Jtheta <- lapply(transformed_par[[2]], function(x) pacf2ArWithJacobian(x[-1], TRUE)$J)
Jtheta <- lapply(transformed_par[[2]], function(x) pacf2ArWithJacobian(x[-1], TRUE)$J)
Judelta <- lapply(transformed_par[[3]], function(x) pacf2ArWithJacobian(x[-1])$J )
JphiJtheta <- dbind(Jphi, Jtheta)
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - There are fewer variables in symmetric approach, so J must be smaller
if(use.symm){
uParLen <- length(unlist(par.skeleton[["udelta"]]))
uPar_seq <- seq_len(uParLen)
Judelta <- lapply(Judelta, function(x) x[uPar_seq, uPar_seq])
JphiJthetaJudelta <- dbind(Jphi, Jtheta, Judelta)
fixed_udelta_symm <- lapply(e_udelta, function(x) x$fixed_symm)
fixed_coef_symm <- unlist(c(fixed_phi, fixed_theta, fixed_udelta_symm,
Par_treg_fixed, Par_const_fixed))
if(length(fixed_coef_symm) == 0)
fixed_coef_symm <- logical(0)
J <- dbind(JphiJthetaJudelta, diag(length(Par_treg)), diag(length(Par_const)))
if(any(fixed_coef_symm))
J <- J[!fixed_coef_symm, !fixed_coef_symm]
}else{
JphiJthetaJudelta <- dbind(Jphi, Jtheta, Judelta)
# TODO: for ss.method = "sarima" this is done by calling again optim().
# see the comments there.
# if(tanh.flag){
# ## tanh applied before sending params to optim()
# ## note that some methods process this internally, e.g. "base";
# ## the correction here is needed only for methods who explicitly request it.
# ## see also fitArma0Model()
# # Jtanh <- diag(1/cosh(arma0_fit$par[seq_len(p+q)])^2)
# # J <- JphiJtheta %*% Jtanh
# phithetapar_only <- c(phipar_only, thetapar_only)
# JphiJtheta <- JphiJtheta / rep(cosh(phithetapar_only)^2, each = length(phithetapar_only))
# stop("unfinished feature")
# }
# J <- dbind(Jphi, Jtheta, diag(length(Par_treg)), diag(length(Par_const)))
# J <- dbind(JphiJtheta, diag(length(Par_treg)), diag(length(Par_const)))
J <- dbind(JphiJthetaJudelta, diag(length(Par_treg)), diag(length(Par_const)))
if(any(fixed_coef))
J <- J[!fixed_coef, !fixed_coef]
}
var_coef <- if(length(coef_hessian) == 0)
matrix(0, 0, 0)
else if(lik.method %in% c("fkf", "kfas"))
# fkf returns the hessian scaled, J H^(-1) J'
# (no, it is from optim! Is the likelihood scaled differently?)
# see the examples in FKF
hessian2vcov(coef_hessian, 1, J)
else
hessian2vcov(coef_hessian, n.used, J)
# This code adds rows and columns (of zeroes) for the fixed params.
# Removing since this is only trouble (and e.g. arima() doesn't do that).
# if(any(fixed_coef)){
# vcov <- matrix(0, length(coef), length(coef))
# vcov[!fixed_coef, !fixed_coef] <- var_coef
# }else
# vcov <- var_coef
vcov <- var_coef
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - Expand out the variance matrix for the repeated parameters
if(use.symm){
.expand_vcov_symm <- function(vcov, par.skeleton){
## where are each of the parameters?
phiInd <- grepl("phi", names(unlist(par.skeleton)))
thetaInd <- grepl("theta", names(unlist(par.skeleton)))
uInd <- grepl("udelta", names(unlist(par.skeleton)))
xregInd <- grepl("xreg", names(unlist(par.skeleton)))
regxInd <- grepl("regx", names(unlist(par.skeleton)))
## Not interested in stationary or regression distinctions
statInd <- phiInd | thetaInd
regInd <- xregInd | regxInd
## save each block
statmat <- vcov[statInd, statInd]
statUmat <- vcov[statInd, uInd]
statregmat <- vcov[statInd, regInd]
Umat <- vcov[uInd, uInd]
Uregmat <- vcov[uInd, regInd]
regmat <- vcov[regInd, regInd]
## expand the symmetric U parameters
uParLen <- length(unlist(par.skeleton[["udelta"]]))
rev_seq <- rev(seq_len(uParLen))
## if there's an even number of parameters, the middle one isn't repeated
if(Udelta_length %% 2 == 0){
rev_seq <- rev_seq[-1]
}
## build the matrix
topBit <- cbind(statmat, statUmat, statUmat[,rev_seq], statregmat)
upperMiddleBit <- cbind(t(statUmat), Umat, Umat[,rev_seq], Uregmat)
lowerMiddleBit <- upperMiddleBit[rev_seq,]
bottomBit <- t(rbind(statregmat, Uregmat, Uregmat[rev_seq,], regmat))
## and return it all
rbind(topBit, upperMiddleBit, lowerMiddleBit, bottomBit)
}
vcov <- .expand_vcov_symm(vcov, par.skeleton)
}
if(is.null(colnames(vcov))){
rownames(vcov) <- colnames(vcov) <- names(coef)[!fixed_coef]
}
## TODO: complete
compact_arima <- as.integer(c(p = length(Phi), q = length(Theta),
s = 1, P = 0, Q = 0))
#browser()
res <- structure(list(
## fields for class Arima (e.g. of the results of arima())
coef = coef,
sigma2 = sigma2,
var.coef = vcov, # var,
mask = mask,
loglik = value_loglik,
aic = aic,
arma = compact_arima, # arma,
residuals = resid,
call = match.call(),
series = y,
code = res_opt$convergence,
n.cond = ncond,
nobs = n.used,
model = ss.model,
## additional fields from me
res_opt = res_opt,
internal = list(
phi = phi,
theta = theta,
delta = delta,
phi_poly = phi_poly,
theta_poly = theta_poly,
delta_poly = delta_poly,
Phi = Phi,
Theta = Theta,
Delta = Delta,
fixed = all.fixed,
nonfixed = coef.nonfixed,
flat_par = flat_par,
fromKalman = fromKalman,
fromKalmanRun = if(exists("fromKalmanRun")) fromKalmanRun else NULL,
transformed_par = transformed_par,
n.coef = n.coef,
trendMaker = trmake,
regxMaker = regxmake
)
## for now inherit from class Arima (the class from arima()),
## eventually do not do this
), class = c("Sarima", "Arima"))
#browser()
res
}
GeoBosh/sarima documentation built on March 27, 2024, 6:31 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions sarimat sarima
Documented in sarima
## Do not edit this file manually.
## It has been automatically generated from *.org sources.
sarima <- function(model, data = NULL, ss.method = "sarima", use.symmetry = FALSE,
SSinit = "Rossignol2011"){
# TODO: for now lik.method is simply "ss.method"
# Need to think over the user-facing arguments for methods.
## N
lik.method <- ss.method
if(!inherits(model, "formula"))
stop("unsupported class ", class(model), " of argument 'model'")
toPolyCoef <- function(item){
switch(item$name,
ar = { phi <<- c( phi , list(item$coef) ) },
ma = { theta <<- c( theta , list(item$coef) ) },
sar = { phi <<- c( phi , list(item$coef) ) },
sma = { theta <<- c( theta , list(item$coef) ) },
i = { delta <<- c( delta , list(item$coef) ) },
si = { delta <<- c( delta , list(item$coef) ) },
s = { delta <<- c( delta , list(item$coef) ) },
ss = { delta <<- c( delta , list(item$coef) ) },
# here item$coef is itself a list!
su = { delta <<- c( delta , item$coef ) },
u = { delta <<- c( delta , item$coef ) },
# unit root for estimation:
uar = { udelta <<- c( udelta, list(item$coef) ) },
stop("not supported yet!")
)
}
Fo <- Formula(model)
yname <- all.vars(formula(Fo, rhs = FALSE)) # name of the response variable
if(length(yname) != 1)
stop("currently there should be exactly one response variable")
data.env <- environment(Fo)
if(!is.null(data)){
if(is.list(data)){ # list (including data.frame)
data.env <- list2env(data, parent = data.env)
environment(Fo) <- data.env
# put the response in y.
# TODO: I had forgotten that I do this! Is it a good idea?
data <- model.part(Fo, data = data, lhs = 1)
data <- as.matrix(data) # TODO: should this be as.data.frame?
}else if(is.environment(data))
stop("'data' from class 'environment' not supported")
else{ # data is the time series # TODO: consider multivariate
data.env <- new.env(parent = data.env)
environment(Fo) <- data.env
data.env[[yname]] <- data
mf0 <- model.frame(Fo, lhs = 1, rhs = FALSE)
data <- model.part(Fo, data = mf0, lhs = 1)
data <- as.matrix(data) # TODO: should this be as.data.frame?
}
}else{
# !!! TODO: do this for the above as well !!!
#
# 15/01/2018 - na.action = NULL, otherwise NA's are dropped
# mf0 <- model.frame(Fo, lhs = 1, rhs = FALSE)
mf0 <- model.frame(Fo, lhs = 1, rhs = FALSE, na.action = NULL)
data <- model.part(Fo, data = mf0, lhs = 1)
data <- as.matrix(data) # TODO: should this be as.data.frame?
}
stopifnot(identical(environment(Fo), data.env))
# get nobs - there should be a better way, since
# calling model frame requires that the
# correct response variable is visible. I
# work only with the lhs, so the other
# variables are not needed at this stage
# (as it have to - cs, ets., are set
# using nobs further below).
Fo.0 <- formula(Fo, lhs = 1, rhs = FALSE)
## 15/01/2018 - bug fix: this drops NA's!
## nobs <- nrow(model.frame(Fo.0)) # TODO: model.frame() needs the environment to be set up!
# any other wy to find nobs?
# DONE: Moving stuff seting the environment before this line
# but leaving this note here in case I forget.
nobs <- nrow(model.frame(Fo.0, na.action = NULL))
#browser()
## NOTE: Do not reuse Fo.0 - its environment gets out of sync with Fo.0o after the following.
## (but the data are already there, so maybe no problem.
# insert an environment for the trend stuff
data.env <- new.env(hash = FALSE, parent = environment(Fo))
data.env$t <- seq_len(nobs)
Fo.trend <- formula(Fo, rhs = 1)
Fo.trend.rhs <- formula(Fo, lhs = FALSE, rhs = 1)
trmake <- NULL
if(!is.empty.model(Fo.trend))
# lhs is needed if there is only intercept, otherwise nobs is unknown
# trmake <- trendMaker(Fo.trend.rhs, time = 1:nobs)
trmake <- trendMaker(Fo.trend, time = 1:nobs)
Fo.sarima <- formula(Fo, rhs = 2) # TODO: check for empty
udelta <- delta <- theta <- phi <- list() # initialise variables, will be updated by toPolyCoef()
molist <- SARIMA(Fo.sarima)
lapply(molist, toPolyCoef)
## renaming Xreg to regx
regxmake <- NULL
if(length(Fo)[2] > 2){
Fo.regx <- formula(Fo, rhs = 3)
if(!is.empty.model(Fo.regx))
regxmake <- trendMaker(Fo.regx, time = 1:nobs)
}else{
Fo.regx <- NULL # for clarity; currently not used
}
#browser()
res <- sarimat(data, phi, theta, delta, udelta, trmake = trmake, regxmake = regxmake,
lik.method = lik.method, use.symm = use.symmetry, SSinit = SSinit)
res$call <- match.call() # replace the call to sarimat() with the call to sarima()
res$formula <- model
res$internal$Fo.xreg <- Fo.trend
res$internal$Fo.sarima <- formula(Fo, lhs = FALSE, rhs = 2) # Fo.sarima
res$internal$Fo.regx <- Fo.regx
#browser()
res
}
sarimat <- function(y, phi, theta, delta, udelta, trmake = NULL, regxmake = NULL, lik.method,
use.symm = FALSE, SSinit = c("Rossignol2011", "gnb", "Gardner1980")){
## 2018-10-12
## SSinit_Q0 <- "Rossignol2011"
## SSinit_Q0 <- "Gardner1980"
SSinit_Q0 <- match.arg(SSinit)
makeArimaInternal <- if(SSinit_Q0 == "gnb")
makeArimaGnb
else
makeARIMA
###### JAMIE: changed 23/08/2018 to include a symmetric approach to U(z)
update_params <- function(par, use.symm = FALSE){
## flat_par contains the parameters as seen by optim,
## except that optim only gets the nonfixed parameters
flat_par[nonfixed] <<- par
## the calculation below depend on the "real" parameters,
## so we tanh transform if necessary to get parcor's, then transform the parcor's
## to coefficients.
## this assumes that the sarima parameters are before any others
## (which is the case, but this is are minder not to change it!)
## using integer indices rather than logical ones to avoid extending
## flags.atanh_sarima to the full length of flat_par.
##
## TODO: In some cases this transformation may make sense for other parameters too.
## But let's sort the sarima parameters first
if(tanh.flag){ ## tanh() parameters to get parcor's
wrk <- flat_par
## fixed parameters declared to need transform, will be transformed to and from.
## TODO: fix this
wrk[ind.atanh_sarima] <- tanh(flat_par[ind.atanh_sarima])
par <- wrk[nonfixed]
}
newpars <- relist(par, par.skeleton)
mapply(function(x, newpar){
x$parcor[x$dyn.index] <- newpar
x$a[-1] <- - pacf2Ar(x$parcor[-1]) # 2022-02-14 was: FitAR::PacfToAR(x$parcor[-1])
}, e_phi, newpars[["phi"]], SIMPLIFY = FALSE)
mapply(function(x, newpar){
x$parcor[x$dyn.index] <- newpar
x$a[-1] <- pacf2Ar(x$parcor[-1]) # 2022-02-14 was: FitAR::PacfToAR(x$parcor[-1])
}, e_theta, newpars[["theta"]], SIMPLIFY = FALSE)
if(use.symm){
mapply(function(x, newpar){
x$parcor_symm[x$dyn.index_symm] <- newpar
x$parcor[x$polyInd] <- x$parcor_symm
x$a[-1] <- - pacf2Ar(x$parcor[-1]) # 2022-02-14 was: FitAR::PacfToAR(x$parcor[-1])
x$coef_symm <- x$a[x$paramInd]
}, e_udelta, newpars[["udelta"]], SIMPLIFY = FALSE)
}else{
mapply(function(x, newpar){
x$parcor[x$dyn.index] <- newpar
x$a[-1] <- - pacf2Ar(x$parcor[-1]) # 2022-02-14 was: FitAR::PacfToAR(x$parcor[-1])
}, e_udelta, newpars[["udelta"]], SIMPLIFY = FALSE)
}
Par_treg <<- if(treg.flag) unlist(newpars[["xreg"]])
else numeric(0)
Par_const <<- if(X.flag) unlist(newpars[["regx"]])
else numeric(0)
NULL
}
makePhiTheta <- function(){
Phi <<- - coef(sprod(phi_poly))[-1]
if(length(Phi) != Phi_length)
Phi <<- c(Phi, numeric(Phi_length - length(Phi)))
Theta <<- coef(sprod(theta_poly))[-1]
if(length(Theta) != Theta_length)
Theta <<- c(Theta, numeric(Theta_length - length(Theta)))
Phi <<- - coef(sprod(phi_poly))[-1]
Delta_prod <<- delta_prod * sprod(udelta_poly)
Delta <<- - coef(Delta_prod)[-1]
if(length(Delta) != Delta_length)
Delta <<- c(Delta, numeric(Delta_length - length(Delta)))
}
makeArimaConst <- function(const, ...){
## see Durbin/Koopman, p. 61
## only const for now, since otherwise Z would be time dependent.
mod <- makeArimaInternal(...)
mod$Z <- c(1, mod$Z)
mod$a <- c(const, mod$a)
## not Matrix::bdiag() for a reason
## - package currently doesn't import Matrix
## - needs to be able to handle e.g. a 0-row matrix
mod$P <- rbind(c(const, numeric(ncol(mod$P))), cbind(0, mod$P))
mod$Pn <- rbind(c(0, numeric(ncol(mod$Pn))), cbind(0, mod$Pn))
mod$T <- rbind(c(1, numeric(ncol(mod$T))), cbind(0, mod$T))
## V = RQR', diag_bind(0, V), i.e. non-random intercept
mod$V <- rbind(c(0, numeric(ncol(mod$V))),
cbind(0, mod$V))
mod
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
loglik <- function(par, use.symm = FALSE){
update_params(par, use.symm)
makePhiTheta() # expand polynomials
if(X.flag){
model_from_makeARIMA <<- makeArimaConst(Par_const, Phi, Theta, Delta, SSinit = SSinit_Q0)
}else
model_from_makeARIMA <<- makeArimaInternal(Phi, Theta, Delta, SSinit = SSinit_Q0)
if(treg.flag)
yc <<- y - tregmat %*% Par_treg
#browser()
fromKalman <<- KalmanLike(yc, model_from_makeARIMA)
fromKalman$Lik
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
loglik_sarima <- function(par, use.symm = FALSE){
update_params(par, use.symm)
makePhiTheta() # expand polynomials
if(X.flag){
model_from_makeARIMA <<- makeArimaConst(Par_const, Phi, Theta, Delta, SSinit = SSinit_Q0)
}else
model_from_makeARIMA <<- makeArimaInternal(Phi, Theta, Delta, SSinit = SSinit_Q0)
if(treg.flag)
yc <<- y - tregmat %*% Par_treg
# fromKalman <<- KalmanLike(yc, model_from_makeARIMA)
fromKalman <<- sarima_KalmanLike(yc, model_from_makeARIMA)
#cat("Phi = ", Phi, " Theta = ", Theta, "\n")
#browser()
fromKalman$Lik
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
ss_sarima <- function(par, use.symm = FALSE){
update_params(par, use.symm)
makePhiTheta() # expand polynomials
if(X.flag){
model_from_makeARIMA <<- makeArimaConst(Par_const, Phi, Theta, Delta, SSinit = SSinit_Q0)
}else
model_from_makeARIMA <<- makeArimaInternal(Phi, Theta, Delta, SSinit = SSinit_Q0)
if(treg.flag)
yc <<- y - tregmat %*% Par_treg
# fromKalman <<- KalmanLike(yc, model_from_makeARIMA)
fromKalman <<- sarima_KalmanLike(yc, model_from_makeARIMA)
fromKalman$css
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
loglik_fkf <- function(par, use.symm = FALSE){
npar.fkf <- length(par)
## @jamie TODO: You set it to this:
## update_params(par, use.symm)
## However, I leave the original, since it doesn't pass all parameters:
## TODO: check if this is oversight.
##
## @georgi (10/10/18) This is as it should be.
##
update_params(par[- npar.fkf])
sigma <<- exp(par[npar.fkf])
makePhiTheta() # expand polynomials
## TODO: check X.flag, etc.
ct <- t(tregmat %*% Par_treg)
dt <- t(regxmat %*% Par_const)
fkf.model <<- xarimaxss(ar = Phi, ma = Theta, sigma = sigma,
xreg = ct, regx = dt, delta = Delta)
fromKalman <<- FKF::fkf(a0 = fkf.model$a0, P0 = fkf.model$P0,
dt = fkf.model$dt, ct = fkf.model$ct,
Tt = fkf.model$Tt, Zt = fkf.model$Zt,
HHt = fkf.model$HHt, GGt = fkf.model$GGt,
yt = matrix(y, nrow = 1))
- fromKalman$logLik # TODO: - $logLik, while in "base" it is $Lik !!!
}
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
loglik_kfas <- function(par, use.symm = FALSE){
npar.fkf <- length(par) # last param is log(sigma2) (In FKF it is log(sigma))
## @jamie TODO: leaving as is,see note above
## @georgi (10/10/18): This is as it should be - as above.
update_params(par[- npar.fkf])
# ## TODO: correct the hessian later
# beta <- tanh(par[- npar.fkf])
# cat("beta = ", beta, "\n")
#
# update_params_kfas(beta)
## TODO: is it better to use sigma as parameter for optim()?
## (better scaling?)
sigma2 <<- exp(par[npar.fkf]) # TODO: check - is this sigma2 or sigma squared?
makePhiTheta() # expand polynomials
# if(X.flag){
# model_from_makeARIMA <<- makeArimaConst(Par_const, Phi, Theta, Delta)
# }else
# model_from_makeARIMA <<- makeARIMA(Phi, Theta, Delta)
# if(treg.flag)
# yc <<- y - tregmat %*% Par_treg
#
# fromKalman <<- KalmanLike(yc, model_from_makeARIMA)
#
# fromKalman$Lik
## TODO: integrated, trends etc,
# ct <- t(tregmat %*% Par_treg)
# dt <- t(regxmat %*% Par_const)
# fkf.model <<- xarimaxss(ar = Phi, ma = Theta, sigma = sigma,
# xreg = ct, regx = dt, delta = Delta)
# from makeUpdateFun_arma0() and KFAS examples
# TODO: only ARMA w.o. intercept for now.
wrk <- try(KFAS::SSMarima(Phi, Theta, Q = sigma2), silent = TRUE)
if(inherits(wrk, "try-error")){
# cat("1e100 !!\n")
return(1e100)
}
kfas_model["T", "arima"] <<- wrk$T
kfas_model["R", "arima"] <<- wrk$R
kfas_model["P1", "arima"] <<- wrk$P1
kfas_model["Q", "arima"] <<- wrk$Q
fromKalman <<- list(logLik = logLik(kfas_model))
- fromKalman$logLik # TODO: - $logLik, while in "base" it is $Lik !!!
}
poly_param <- function(x, fixed = attr(x, "fixed"), nseasons = attr(x, "nseasons"),
d = attr(x, "d", exact = TRUE),
atanh.tr = attr(x, "atanh.tr", exact = TRUE) ){
coef <- if(!is.null(attr(x, "sign")) && attr(x, "sign") == "-")
c(1, - x)
else
c(1, x)
# polynom() drops trailing zeroes and doesn't accept NA's.
# to protect against these, create a polynomial of 1's first
res <- polynom(rep(1, length(coef)))
## TODO: krapka - replace NA's with zeroes
if(anyNA(coef))
coef[is.na(coef)] <- 0
e <- environment(res)
e$a <- coef
e$fixed <- if(is.null(fixed))
c(TRUE, rep(FALSE, length(coef) - 1))
else
c(TRUE, fixed)
e$dyn.index <- which(!e$fixed)
## 18/01/2017 - support for atanh transform
## TODO: maybe different defaults for AR and MA? but here it is not known if
## the parameter is AR or MA.
## It is best for the default here to be FALSE since otherwise the caller
## will have to take care of this even for model elements for which this
## doesn't make sense. In any case, the default here is not that important
## since sarima() may set these.
e$atanh.tr <- if(is.null(atanh.tr))
c(FALSE, rep(FALSE, length(coef) - 1))
else
c(FALSE, atanh.tr)
## TODO: for now there seems no need for analog of dyn.index (see above)
if(!is.null(nseasons)){
## TODO: check for whole number
if(is.numeric(nseasons) && length(nseasons) == 1 && nseasons > 0){
e$nseasons <- nseasons
}else
stop("'nseasons' must be a positive integer or NULL")
}
if(!is.null(d) && d > 1){
e$d <- d
}
if(!is.null(dnam <- attr(x, "dispname")))
e$dispname <- dnam
res
}
sprod <- function(pp){ # pp is short for parameter polynomial
f <- function(poly){
s <- environment(poly)$nseasons
if(is.null(s))
poly
else
poly(poly_x^s)
}
prod(as.polylist(lapply(pp, f)))
}
u_sprod <- function(pp){ # pp is short for parameter polynomial
f <- function(poly){
s <- environment(poly)$nseasons
if(!is.null(s))
poly <- poly(poly_x^s)
d <- environment(poly)$d
if(is.null(d)) poly else poly^d
}
prod(as.polylist(lapply(pp, f)))
}
sdegree <- function(pp){
e <- environment(pp)
s <- e$nseasons
degree <- length(e$a) - 1
if(!is.null(s))
degree <- degree * s
if(!is.null(e$d))
degree <- degree * e$d
degree
}
## TODO: this is krapka for naming multple ar() terms, etc.
poly_param_plus <- function(parlist){
no <- 0
sno <- 0
res <- vector(length(parlist), mode = "list")
for(i in seq_along(parlist)){
wrk <- poly_param(parlist[[i]])
e <- environment(wrk)
if(is.null(e$nseasons)){
e$name.counter <- no
no <- no + 1
}else{
e$name.counter <- sno
sno <- sno + 1
}
res[[i]] <- wrk
}
as.polylist(res)
}
# phi, theta and delta are used only here.
# TODO: maybe the they can come as polynomials to sarimat() ?
phi_poly <- poly_param_plus(phi)
theta_poly <- poly_param_plus(theta)
delta_poly <- poly_param_plus(delta)
udelta_poly <- poly_param_plus(udelta)
# Delta is fixed, so needs to be computed only once
## TODO: (1-x)^d is kept as 1-x and a 'd' is a sepaarate varible.
## the element of delta_poly corresponding to this will print as 1 - x.
## u_sprod(), print.Sarima(), summary.Sarima() know this but ing eneral this is very fragile.
delta_prod <- u_sprod(delta_poly)
udelta_prod <- sprod(udelta_poly) # NOTE: this is not u_sprod(udelta_poly)
Par_const <- Par_treg <- Phi <- Theta <- numeric(0)
yc <- y
sigma <- sd(y) # for fkf
sigma2 <- var(y) # for kfas
kfas_model <- fkf.model <- fromKalman <- model_from_makeARIMA <- NULL
Delta_prod <- delta_prod * udelta_prod
Delta <- - coef(Delta_prod)[-1]
e_phi <- lapply(phi_poly, environment)
e_theta <- lapply(theta_poly, environment)
e_delta <- lapply(delta_poly, environment)
e_udelta <- lapply(udelta_poly, environment)
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - add extra arguments to the environment to help this approach
if(use.symm){
.nonstationaryEdit <- function(e){
## Parcors will be symmetric inside the extreme parcors
## so we only want half the parameters
len <- length(e$a)
midpoint <- ceiling(len/2 - 1)
mid_seq <- seq_len(midpoint)
## fixed and atanh.tr, accounting for symmetry
coef_a <- e$a[-c(1,len)]
fixed_a <- e$fixed[-c(1,len)]
atanh_a <- e$atanh.tr[-c(1,len)]
fixed_symm <- fixed_a[mid_seq] | rev(fixed_a)[mid_seq]
atanh_symm <- atanh_a[mid_seq] | rev(atanh_a)[mid_seq]
## where to find b parameters
paramInd <- mid_seq + 1
## starting coefficients
coef_symm <- e$a[-c(1,len)][mid_seq]
## if any value is fixed, make sure we take it appropriately
if(any(fixed_symm)){
fix_a <- which(fixed_a)
for(i in seq_along(fix_a)){
if(fix_a[i] > midpoint){
fix_symm <- len - 2 - fix_a[i]
paramInd[fix_symm] <- fix_a[i] + 1
}else{
fix_symm <- fix_a[i]
}
coef_symm[fix_symm] <- coef_a[fix_a[i]]
}
}
## add to environment
e$fixed_symm <- fixed_symm
e$atanh.tr_symm <- atanh_symm
e$dyn.index_symm <- which(!fixed_symm)
e$paramInd <- paramInd
e$coef_symm <- coef_symm
## where b parameters can be found
e$polyInd <- c(seq_along(coef_symm) + 1,
len - seq_along(coef_symm))
NULL
}
lapply(e_udelta, .nonstationaryEdit)
}
e_phi_theta <- c(e_phi, e_theta)
e_phi_theta_delta <- c(e_phi_theta, e_delta)
e_phi_theta_udelta <- c(e_phi_theta, e_udelta)
Phi_length <- if(length(phi_poly) > 0) sum(sapply(phi_poly, sdegree)) else 0
Theta_length <- if(length(theta_poly) > 0) sum(sapply(theta_poly, sdegree)) else 0
Udelta_length <- if(length(udelta_poly) > 0) sum(sapply(udelta_poly, sdegree)) else 0
delta_length <- if(length(delta_poly) > 0) sum(sapply(delta_poly, sdegree)) else 0
Delta_length <- delta_length + Udelta_length
prepare_parcor <- function(envir){
if(is.null(envir$parcor)){
envir$parcor <- c(1, ar2Pacf(- envir$a[-1]))
}else if(length(envir$parcor) != length(envir$a))
stop("lengths of coef and parcor are not equal")
if(anyNA(envir$parcor)) # todo: give a warning in this case?
envir$parcor[is.na(envir$parcor)] <- 0
NULL
}
# set initial values for parcor's (transformed params)
# also: "parcor's don't make sense for delta!
lapply(e_phi_theta, prepare_parcor)
lapply(e_udelta, function(e) e$parcor <- e$a) # TODO: $parcor should be set outright, this is a hack
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - update the parameters that will be driving the equation
if(use.symm) lapply(e_udelta, function(e) e$parcor_symm <- e$parcor[e$paramInd])
# TODO: this is just wrong, but not used:
lapply(e_delta, function(e) e$parcor <- e$a) # TODO:delta !!!!!
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - The symmetric approach uses fewer parameters, so the length
###### of variables in the functions defined at the start needs to
###### be the correct size, etc.
if(use.symm){
flags.atanh_phi_theta <- lapply(e_phi_theta, function(x)
if(isTRUE(x$operator)) logical(0) else x$atanh.tr)
flags.atanh_udelta <- lapply(e_udelta, function(x)
if(isTRUE(x$operator)) logical(0) else x$atanh.tr_symm)
flags.atanh_sarima <- unlist(c(flags.atanh_phi_theta, flags.atanh_udelta))
if(length(flags.atanh_sarima) == 0) ## in case it is list()
flags.atanh_sarima <- logical(0)
ind.atanh_sarima <- which(flags.atanh_sarima)
tanh.flag <- any(flags.atanh_sarima)
Par_phi_theta <- lapply(e_phi_theta, function(x)
if(isTRUE(x$operator)) numeric(0) else x$parcor)
Par_udelta <- lapply(e_udelta, function(x)
if(isTRUE(x$operator)) numeric(0) else x$parcor_symm)
Par_sarima <- unlist(c(Par_phi_theta, Par_udelta))
fixed_phi_theta <- lapply(e_phi_theta, function(x)
if(isTRUE(x$operator)) logical(0) else x$fixed)
fixed_udelta <- lapply(e_udelta, function(x)
if(isTRUE(x$operator)) logical(0) else x$fixed_symm)
fixed <- unlist(c(fixed_phi_theta, fixed_udelta))
}else{
flags.atanh_sarima <- unlist(lapply(e_phi_theta_udelta,
function(x)
if(isTRUE(x$operator)) logical(0) else x$atanh.tr
))
if(length(flags.atanh_sarima) == 0) ## in case it is list()
flags.atanh_sarima <- logical(0)
ind.atanh_sarima <- which(flags.atanh_sarima)
tanh.flag <- any(flags.atanh_sarima)
Par_sarima <- unlist(lapply(e_phi_theta_udelta,
function(x)
if(isTRUE(x$operator)) numeric(0) else x$parcor
))
## @jamie TODO: !!!
## You have changed the logic even n the non-sym case, here and at other places.
## It is extremely dangerous to do multiple unrelated sweeping changes
## at the same time, especcially without proper testing.
## e.g. this assignement has been mmoved from further below here, silently.
fixed <- unlist(lapply(e_phi_theta_udelta,
function(x)
if(isTRUE(x$operator)) logical(0) else x$fixed
))
}
stopifnot(length(Par_sarima) == length(flags.atanh_sarima))
if(tanh.flag)
Par_sarima[flags.atanh_sarima] <- atanh(Par_sarima[flags.atanh_sarima])
X.flag <- !is.null(regxmake)
treg.flag <- !is.null(trmake)
## now prepare the treg parameters
## time regression
## set up xreg
if(treg.flag){
tregmat <- trmake()
if(length(Delta) > 0){
treg.temp <- trmake((1 - length(Delta)) : nrow(tregmat))
treg.temp <- cbind(c(rep(NA_real_, length(Delta)), y), treg.temp)
# TODO: use filter()
## start from the last row to do it in-place
for(i in nrow(treg.temp):(length(Delta) + 1)){
for(j in 1:length(Delta)){
treg.temp[i, ] <- treg.temp[i, ] - Delta[j] * treg.temp[i - j, ]
}
}
## drop the added rows
treg.temp <- treg.temp[-(1:length(Delta)), ]
ydiff <- treg.temp[-(1:length(Delta)), 1]
tregdiff <- treg.temp[-(1:length(Delta)), -1]
fitinit <- lm(ydiff ~ tregdiff - 1, na.action = na.omit)
Par_treg <- coef(fitinit)
}else{
# Tova nyama da e vyarno ako ima regx:
# if(ncol(tregmat) == 1 && all(tregmat == 1)){
# colnames(tregmat) <- "mean"
# }
fitinit <- lm(y ~ tregmat - 1, na.action = na.omit)
Par_treg <- coef(fitinit)
}
}else{#treg.flag = FALSE
tregmat <- matrix(0, nrow = length(y), ncol = 0)
Par_treg <- numeric(0)
}
## set up regx
if(X.flag){
regxmat <- regxmake()
fitinit <- lm(y ~ regxmat - 1, na.action = na.omit)
Par_const <- coef(fitinit)
}else{
regxmat <- matrix(0, nrow = length(y), ncol = 0)
Par_const <- numeric(0)
}
# note: flat_par includes the coeffcients of zero power of the polynomials
flat_par <- c(Par_sarima, Par_treg, Par_const)
# TODO: for now no fixed treg parameters
# use 'offset' if possible?
Par_treg_fixed <- rep(FALSE, length(Par_treg))
Par_const_fixed <- rep(FALSE, length(Par_const))
fixed <- c(fixed, Par_treg_fixed, Par_const_fixed)
## note: flat_par includes the coeffcients of zero power of the polynomials;
## so, flat_par[nonfixed] includes them, too!
nonfixed <- !fixed
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - par.skeleton[["udelta]] will have fewer dynamic variables
if(use.symm){
par.skeleton <- list(phi = lapply(e_phi, function(e) e$dyn.index),
theta = lapply(e_theta, function(e) e$dyn.index),
udelta = lapply(e_udelta, function(e) e$dyn.index_symm),
# delta = lapply(e_delta, function(e) e$dyn.index),
xreg = list(seq_along(Par_treg)),
regx = list(seq_along(Par_const)))
}else{
par.skeleton <- list(phi = lapply(e_phi, function(e) e$dyn.index),
theta = lapply(e_theta, function(e) e$dyn.index),
udelta = lapply(e_udelta, function(e) e$dyn.index),
# delta = lapply(e_delta, function(e) e$dyn.index),
xreg = list(seq_along(Par_treg)),
regx = list(seq_along(Par_const)))
}
## optimise
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - each method now includes an argument to use a symmetric approach
switch(lik.method,
base = {
loglik(flat_par[nonfixed], use.symm = use.symm)
res_opt <- optim(flat_par[nonfixed], loglik, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
# update the quantities with the result from the optim:
# (the last call of loglik() by optim() may not be the final result)
loglik(res_opt$par, use.symm = use.symm) # identical(res_opt$value, fromKalman$Lik) should be TRUE
},
sarima = {
loglik_sarima(flat_par[nonfixed], use.symm = use.symm)
res_opt <- optim(flat_par[nonfixed], loglik_sarima, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
if(tanh.flag){
## lazy: call optim with the untransformed parameters with zero maxit,
## simply to get the hessian for the untransformed parameters.
## todo: however the current setup makes this difficult
## save some vars
tanh.flag <- FALSE
sav_flat_par <- flat_par
## tanh() to get parcor's
flat_par[ind.atanh_sarima] <- tanh(flat_par[ind.atanh_sarima])
## this is only needed for the hessian
res_opt_alt <- optim(flat_par[nonfixed], loglik_sarima, use.symm = use.symm,
method = "BFGS", hessian = TRUE,
control = list(maxit = 0L))
flat_par <- sav_flat_par # restore saved vars
tanh.flag <- TRUE
}
# update the quantities with the result from the optim:
# (the last call of loglik() by optim() may not be the final result)
loglik_sarima(res_opt$par, use.symm = use.symm) # identical(res_opt$value, fromKalman$Lik) should be TRUE
},
css = {
ss_sarima(flat_par[nonfixed], use.symm = use.symm)
res_opt <- optim(flat_par[nonfixed], ss_sarima, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
if(tanh.flag){
## save some vars
tanh.flag <- FALSE
sav_flat_par <- flat_par
## tanh() to get parcor's
flat_par[ind.atanh_sarima] <- tanh(flat_par[ind.atanh_sarima])
## this is only needed for the hessian
res_opt_alt <- optim(flat_par[nonfixed], ss_sarima, use.symm = use.symm,
method = "BFGS", hessian = TRUE,
control = list(maxit = 0L))
flat_par <- sav_flat_par # restore saved vars
tanh.flag <- TRUE
}
ss_sarima(res_opt$par, use.symm = use.symm)
## 10/10/18 method changed to 'sarima' to avoid duplicating code later
lik.method <- "sarima"
},
fkf = {
loglik_fkf(c(flat_par[nonfixed], log(sigma)), use.symm = use.symm)
res_opt <- optim(c(flat_par[nonfixed], log(sigma)),
loglik_fkf, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
loglik_fkf(res_opt$par, use.symm = use.symm)
},
kfas = {
# KFAS::SSMarima() returns a list
# kfas_model <- try(KFAS::SSMarima(Phi, Theta, Q = sigma2), silent = TRUE)
## loglik_kfas() doesn't recreate the full model,
## so create it in advance (could add if(is.null(kfas_model)) ... to loglik_kfas)
update_params(flat_par[nonfixed], use.symm = use.symm)
makePhiTheta() # expand polynomials
## 2022-01-19 moved KFAS to Suggests.
##
## however, then SSMarima is not visible and in a model formula can't
## be qualified with 'KFAS::'. So, copying it here.
## TODO: move this assignment elsewhere (here it is called every
## time the likelihood is computed, the cost is negligible though)
SSMarima <- KFAS::SSMarima
kfas_model <- KFAS::SSModel(y ~ -1 + SSMarima(Phi, Theta, Q = sigma2), H = 0)
# kfas_init <- c(atanh(flat_par[nonfixed]), log(sigma2))
kfas_init <- c(flat_par[nonfixed], log(sigma2))
loglik_kfas(kfas_init, use.symm = use.symm)
res_opt <- optim(kfas_init, loglik_kfas, use.symm = use.symm,
method = "BFGS", hessian = TRUE)
loglik_kfas(res_opt$par, use.symm = use.symm)
},
arima = {
## this methods calls arima()
if(length(Par_const) > 0)
stop("method 'arima' can handle 'xreg', but not 'regx'")
## TODO: finish this
stop("method arima is unfinished")
},
stop("invalid 'lik.method'")
)
fcoef <- function(x, minus = FALSE){
res <- x$a[-1]
if(minus)
res <- - res
if(!is.null(nam <- x$dispname)){
suffix <- letters[x$name.counter] # should be ok for NULL and 0
names(res) <- paste0(nam, seq_along(res), suffix)
}
res
}
coef_phi <- lapply(e_phi, function(x) fcoef(x, TRUE))
fixed_phi <- lapply(e_phi, function(x) x$fixed[-1])
coef_theta <- lapply(e_theta, fcoef)
fixed_theta <- lapply(e_theta, function(x) x$fixed[-1])
coef_delta <- lapply(e_delta, function(x) fcoef(x, TRUE))
fixed_delta <- lapply(e_delta, function(x) x$fixed[-1])
coef_udelta <- lapply(e_udelta, function(x) fcoef(x, TRUE))
fixed_udelta <- lapply(e_udelta, function(x) x$fixed[-1])
all.coef <- unlist(c(coef_phi, coef_theta, coef_udelta, coef_delta, Par_treg, Par_const))
all.fixed <- unlist(c(fixed_phi, fixed_theta, fixed_udelta, fixed_delta, Par_treg_fixed, Par_const_fixed))
coef.nonfixed <- unlist(all.coef)[!all.fixed]
# TODO: think about better names, can clean them up here
if(length(Par_treg) > 0)
names(Par_treg) <- colnames(tregmat)
if(length(Par_const) > 0)
names(Par_const) <- colnames(regxmat)
# no delta here, but all this needs further thought
# TODO: need to add numbers to c("ar", "ar"), etc
coef <- unlist(c(coef_phi, coef_theta, coef_udelta, Par_treg, Par_const))
fixed_coef <- unlist(c(fixed_phi, fixed_theta, fixed_udelta, Par_treg_fixed, Par_const_fixed))
## TODO: Tova e krapka, ako e list() ili NULL , !fixed_coef po-dolu dava greshka.
if(length(fixed_coef) == 0)
fixed_coef <- logical(0)
n.coef <- c(arma = length(unlist(c(coef_phi, coef_theta))),
# TODO: reinstate but change predict.Sarima to use indexing with names
uar = length(unlist(coef_udelta)),
xreg = length(Par_treg),
regx = length(Par_const))
mask <- !fixed_coef
## n.cond in arima is only for method css, maybe should drop it here.
ncond <- 0 # max(length(Phi), length(Theta) + 1) # TODO: Fix this!
# in arima() n.used is:
# n.used <- sum(!is.na(x)) - length(Delta)
# (with additional correction for missing values in xreg).
# TODO: do the correction for xreg.
notna <- !is.na(y)
n.used <- if(length(Delta) == 0)
sum(notna)
else if(all(!is.na(y[1:length(Delta)])))
sum(notna) - length(Delta)
else{ # NA's among the first length(Delta) values - additional loss
# if there is a stretch of length(Delta) non-NA's all subsequent obs. are ok
# but the exact number of non-ok depends also on zero coefficients in Delta
# (e.g. for 1 - B^12)
## TODO: remove this or further modify at least for ss.method = "sarima"
available <- !is.na(y)
statable <- rep(FALSE, length(y))
ddind <- which(Delta != 0)
okseq <- 0
for(t in (length(Delta) + 1):length(y)){
statable[t] <- all(statable[t - ddind] | available[t - ddind])
if(statable[t])
okseq <- okseq + 1
else
okseq <- 0
if(okseq >= length(Delta)){
if(t < length(y))
statable[(t + 1):length(y)] <- TRUE
break
}
}
## !!! TODO: statable may be useful to choose which residuals to use etc.
sum(statable)
}
switch(lik.method,
base = {
# arima(): value <- 2 * n.used * res$value + n.used + n.used * log(2 * pi)
# aic <- if(method != "css") value + 2*sum(mask) + 2 else NA
# KalmanLike() returns minus log-likelihood, upto a scaling and addition of constants
# (not fully documented, the following is a guess from trial and error and regression)
# equivalently, can replace fromKalman$Lik with res_opt$value
value_loglik <- - n.used * fromKalman$Lik - 0.5 * n.used * log(2 * pi) - 0.5 * n.used
# I don't fully understand argument 'nit', but setting it to zero seems
# to be right for prediction (setting it to -1 gives different results), compare:
# fromKalmanRun <- KalmanRun(y, model_from_makeARIMA, -1, update = TRUE)
# fromKalmanRun <- KalmanRun(y, model_from_makeARIMA, 0, update = TRUE)
fromKalmanRun <- KalmanRun(yc, model_from_makeARIMA, 0, update = TRUE)
resid <- fromKalmanRun$resid
### from jamie: these are standardised residuals!
# TODO: sigma2 needs more care.
# This reproduces (it seems) arima()'s sigma2 in the absence of NA's
ndrop_sigma2 <- length(Delta)
sigma2 <- if(ndrop_sigma2 == 0)
mean(resid^2, na.rm = TRUE)
else
mean(resid[-seq_len(ndrop_sigma2)]^2, na.rm = TRUE)
sigma <- sqrt(sigma2)
## should be ready for prediction (minus exogeneous vars)
ss.model <- attr(fromKalmanRun, "mod") # model_from_makeARIMA
coef_hessian <- res_opt$hessian
#browser()
},
sarima = {
## for now same as "base" but calling my likelihood code
value_loglik <- - n.used * fromKalman$Lik - 0.5 * n.used * log(2 * pi) - 0.5 * n.used
# I don't fully understand argument 'nit', but setting it to zero seems
# to be right for prediction (setting it to -1 gives different results), compare:
# fromKalmanRun <- KalmanRun(y, model_from_makeARIMA, -1, update = TRUE)
# fromKalmanRun <- KalmanRun(y, model_from_makeARIMA, 0, update = TRUE)
## TODO: this needs changing, but lets first sort out the estimation
## The C++ function uniKalmanLikelihood0b needs amending to return states and
## the updated model
## fromKalmanRun <- sarima_KalmanRun(yc, model_from_makeARIMA, 0, update = TRUE)
## (2018-08-23) TODO:
## @jamie: In such cases please put a note has happened here.
## Is this a follow up to the aboe note, if so to what extent, etc.
## I assume that you have tested the equivalence of these two in the context.
## 2018-08-23 was:
## fromKalmanRun <- KalmanRun(yc, model_from_makeARIMA, 0, update = TRUE)
fromKalmanRun <- sarima_KalmanRun(yc, model_from_makeARIMA, 0, update = TRUE)
resid <- fromKalmanRun$resid
# TODO: sigma2 needs more care.
# This reproduces (it seems) arima()'s sigma2 in the absence of NA's
ndrop_sigma2 <- length(Delta)
sigma2 <- if(ndrop_sigma2 == 0)
mean(resid^2, na.rm = TRUE)
else
mean(resid[-seq_len(ndrop_sigma2)]^2, na.rm = TRUE)
sigma <- sqrt(sigma2)
## should be ready for prediction (minus exogeneous vars)
ss.model <- attr(fromKalmanRun, "mod") # model_from_makeARIMA
## TODO: it may be better not to call again optim
## but to transform the hessian as in fitArma0Model()
coef_hessian <- if(tanh.flag)
res_opt_alt$hessian
else
res_opt$hessian
},
fkf = {
value_loglik <- fromKalman$logLik
## fkf puts other useful stuff in the answers
resid <- as.vector(fromKalman$vt)
sigma2 <- sigma^2
ss.model <- fkf.model
nrhess <- dim(res_opt$hessian)[1] # hess is square matrix
coef_hessian <- res_opt$hessian[- nrhess, - nrhess] # drop logsigma
},
kfas = {
value_loglik <- fromKalman$logLik
## fkf puts other useful stuff in the answers
resid <- as.vector(fromKalman$vt) # TODO: this is incorrect!
ss.model <- kfas_model
nrhess <- dim(res_opt$hessian)[1] # hess is square matrix
coef_hessian <- res_opt$hessian[- nrhess, - nrhess] # drop logsigma
},
arima = {
stop("Should not get this far, 'arima' method not supported yet.")
},
stop("'lik.method' not found.")
)
aic <- - 2 * value_loglik + 2 * sum(nonfixed) + 2
transformed_par <- list(lapply(e_phi, function(x) x$parcor),
lapply(e_theta, function(x) x$parcor),
lapply(e_udelta, function(x) x$parcor),
lapply(e_delta, function(x) x$parcor) )
phipar_only <- unlist(lapply(e_phi, function(x) x$a[-1]))
phipar_fixed <- unlist(lapply(e_phi, function(x) x$fixed[-1]))
thetapar_only <- unlist(lapply(e_theta, function(x) x$a[-1]))
thetapar_fixed <- unlist(lapply(e_theta, function(x) x$fixed[-1]))
# Jphi, Jtheta are lists of matrices
# transformed par are the parcor's
Jphi <- lapply(transformed_par[[1]], function(x) pacf2ArWithJacobian(x[-1])$J )
# Jtheta <- lapply(transformed_par[[2]], function(x) pacf2ArWithJacobian(x[-1], TRUE)$J)
Jtheta <- lapply(transformed_par[[2]], function(x) pacf2ArWithJacobian(x[-1], TRUE)$J)
Judelta <- lapply(transformed_par[[3]], function(x) pacf2ArWithJacobian(x[-1])$J )
JphiJtheta <- dbind(Jphi, Jtheta)
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - There are fewer variables in symmetric approach, so J must be smaller
if(use.symm){
uParLen <- length(unlist(par.skeleton[["udelta"]]))
uPar_seq <- seq_len(uParLen)
Judelta <- lapply(Judelta, function(x) x[uPar_seq, uPar_seq])
JphiJthetaJudelta <- dbind(Jphi, Jtheta, Judelta)
fixed_udelta_symm <- lapply(e_udelta, function(x) x$fixed_symm)
fixed_coef_symm <- unlist(c(fixed_phi, fixed_theta, fixed_udelta_symm,
Par_treg_fixed, Par_const_fixed))
if(length(fixed_coef_symm) == 0)
fixed_coef_symm <- logical(0)
J <- dbind(JphiJthetaJudelta, diag(length(Par_treg)), diag(length(Par_const)))
if(any(fixed_coef_symm))
J <- J[!fixed_coef_symm, !fixed_coef_symm]
}else{
JphiJthetaJudelta <- dbind(Jphi, Jtheta, Judelta)
# TODO: for ss.method = "sarima" this is done by calling again optim().
# see the comments there.
# if(tanh.flag){
# ## tanh applied before sending params to optim()
# ## note that some methods process this internally, e.g. "base";
# ## the correction here is needed only for methods who explicitly request it.
# ## see also fitArma0Model()
# # Jtanh <- diag(1/cosh(arma0_fit$par[seq_len(p+q)])^2)
# # J <- JphiJtheta %*% Jtanh
# phithetapar_only <- c(phipar_only, thetapar_only)
# JphiJtheta <- JphiJtheta / rep(cosh(phithetapar_only)^2, each = length(phithetapar_only))
# stop("unfinished feature")
# }
# J <- dbind(Jphi, Jtheta, diag(length(Par_treg)), diag(length(Par_const)))
# J <- dbind(JphiJtheta, diag(length(Par_treg)), diag(length(Par_const)))
J <- dbind(JphiJthetaJudelta, diag(length(Par_treg)), diag(length(Par_const)))
if(any(fixed_coef))
J <- J[!fixed_coef, !fixed_coef]
}
var_coef <- if(length(coef_hessian) == 0)
matrix(0, 0, 0)
else if(lik.method %in% c("fkf", "kfas"))
# fkf returns the hessian scaled, J H^(-1) J'
# (no, it is from optim! Is the likelihood scaled differently?)
# see the examples in FKF
hessian2vcov(coef_hessian, 1, J)
else
hessian2vcov(coef_hessian, n.used, J)
# This code adds rows and columns (of zeroes) for the fixed params.
# Removing since this is only trouble (and e.g. arima() doesn't do that).
# if(any(fixed_coef)){
# vcov <- matrix(0, length(coef), length(coef))
# vcov[!fixed_coef, !fixed_coef] <- var_coef
# }else
# vcov <- var_coef
vcov <- var_coef
###### JAMIE (23/08/2018): include a symmetric approach to U(z)
###### - Expand out the variance matrix for the repeated parameters
if(use.symm){
.expand_vcov_symm <- function(vcov, par.skeleton){
## where are each of the parameters?
phiInd <- grepl("phi", names(unlist(par.skeleton)))
thetaInd <- grepl("theta", names(unlist(par.skeleton)))
uInd <- grepl("udelta", names(unlist(par.skeleton)))
xregInd <- grepl("xreg", names(unlist(par.skeleton)))
regxInd <- grepl("regx", names(unlist(par.skeleton)))
## Not interested in stationary or regression distinctions
statInd <- phiInd | thetaInd
regInd <- xregInd | regxInd
## save each block
statmat <- vcov[statInd, statInd]
statUmat <- vcov[statInd, uInd]
statregmat <- vcov[statInd, regInd]
Umat <- vcov[uInd, uInd]
Uregmat <- vcov[uInd, regInd]
regmat <- vcov[regInd, regInd]
## expand the symmetric U parameters
uParLen <- length(unlist(par.skeleton[["udelta"]]))
rev_seq <- rev(seq_len(uParLen))
## if there's an even number of parameters, the middle one isn't repeated
if(Udelta_length %% 2 == 0){
rev_seq <- rev_seq[-1]
}
## build the matrix
topBit <- cbind(statmat, statUmat, statUmat[,rev_seq], statregmat)
upperMiddleBit <- cbind(t(statUmat), Umat, Umat[,rev_seq], Uregmat)
lowerMiddleBit <- upperMiddleBit[rev_seq,]
bottomBit <- t(rbind(statregmat, Uregmat, Uregmat[rev_seq,], regmat))
## and return it all
rbind(topBit, upperMiddleBit, lowerMiddleBit, bottomBit)
}
vcov <- .expand_vcov_symm(vcov, par.skeleton)
}
if(is.null(colnames(vcov))){
rownames(vcov) <- colnames(vcov) <- names(coef)[!fixed_coef]
}
## TODO: complete
compact_arima <- as.integer(c(p = length(Phi), q = length(Theta),
s = 1, P = 0, Q = 0))
#browser()
res <- structure(list(
## fields for class Arima (e.g. of the results of arima())
coef = coef,
sigma2 = sigma2,
var.coef = vcov, # var,
mask = mask,
loglik = value_loglik,
aic = aic,
arma = compact_arima, # arma,
residuals = resid,
call = match.call(),
series = y,
code = res_opt$convergence,
n.cond = ncond,
nobs = n.used,
model = ss.model,
## additional fields from me
res_opt = res_opt,
internal = list(
phi = phi,
theta = theta,
delta = delta,
phi_poly = phi_poly,
theta_poly = theta_poly,
delta_poly = delta_poly,
Phi = Phi,
Theta = Theta,
Delta = Delta,
fixed = all.fixed,
nonfixed = coef.nonfixed,
flat_par = flat_par,
fromKalman = fromKalman,
fromKalmanRun = if(exists("fromKalmanRun")) fromKalmanRun else NULL,
transformed_par = transformed_par,
n.coef = n.coef,
trendMaker = trmake,
regxMaker = regxmake
)
## for now inherit from class Arima (the class from arima()),
## eventually do not do this
), class = c("Sarima", "Arima"))
#browser()
res
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.