[Robkalman-commits] r83 - pkg/robKalman/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 15 16:34:21 CEST 2025
Author: ruckdeschel
Date: 2025-05-15 16:34:21 +0200 (Thu, 15 May 2025)
New Revision: 83
Removed:
pkg/robKalman/R/arGM.R
pkg/robKalman/R/arGMinternal.R
Log:
arGM moved to pkg robspec
Deleted: pkg/robKalman/R/arGM.R
===================================================================
--- pkg/robKalman/R/arGM.R 2025-05-15 14:31:23 UTC (rev 82)
+++ pkg/robKalman/R/arGM.R 2025-05-15 14:34:21 UTC (rev 83)
@@ -1,114 +0,0 @@
-arGM <- function (x, order=1,
- chr=1.5, iterh=maxiter, cbr=5.0, iterb=maxiter,
- psi2="Tukey", c=4.0, type="Mallows",
- k=1.5, maxiter=100, tol=1e-08, equal.LS=FALSE, ...)
-{
-##################################
-##
-## R-function: arGM - GM-estimates for AR parameters
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-26)
-## References:
-## [MarZ78] R.D. Martin & J.E. Zeh, Generalized M-estimates for
-## Autoregression Including Small-sample Efficiency
-## Robustness (1978)
-## [Mart80] R.D. Martin, Robust Estimation of Autoregressive Models (1980)
-## [MarT82b] R.D. Martin & D.J. Thomson, Robust-resistent Spectrum
-## Estimation (1982)
-## [StoDut87] N. Stockinger & R. Dutter, Robust Time Series Analysis:
-## A Survey (1987)
-##
-##################################
-
-## Parameters:
-## x ... univarite time series (vector)
-## order ... order of AR(p) process (default: order=1)
-## chr ... tuning konstant for Huber's psi function (default: chr=1.5)
-## iterh ... number of iterations for IWLS-alogrithm using
-## Huber's psi function (default: iterh=maxiter)
-## cbr ... tuning konstant for Tukey's psi function (default: cbr=5.0)
-## iterb ... number of iterations for IWLS-alogrithm using
-## Tukey's psi function (default: iterb=maxiter)
-## psi2 ... influence function to determine the 'largness of z_i',
-## either "Ident", "Huber" or "Tukey" (default: "Tukey")
-## c ... tuning constant for psi2 (default: c=4.0)
-## type ... type of GM-estimates, either "Mallows" or "Schweppe"
-## (default: "Mallows")
-## k ... tuning constant for centering (default: k=1.5)
-## maxiter ... maximal number of iteration (default: maxiter=100)
-## tol ... tolerance level (default: tol=1e-08)
-## equal.LS ... logical, for testing purpose only (default: equal.LS=FALSE)
-## ... ... further parameters to be passed to the functions 'HuberM' or
-## 'hubers'
-
- ## Variable definitions:
-
- s <- c()
- Phi <- matrix()
- w <- NA
- BH <- BB <- NA
- niterh <- niterb <- niter.testing <- NA
-
- ## Centering:
-
- ## x.huber <- HuberM(x, ...) # as proposed in [StoDut87]
- x.huber <- hubers(x, ...) # as proposed in [MarZ78], [Mart80]
- x <- x - x.huber$mu
- sx <- x.huber$s
-
- ## Main:
-
- for (p in 1:order) {
- ARmodel <- .ARmodel(x, p)
- y <- ARmodel$y
- Z <- ARmodel$Z
- invCp <- .invCp(p, c(sx, s), Phi)
- Weights <- .Weights(p, Z, invCp, type, psi2, c)
- u <- Weights$u
- v <- Weights$v
- startval <- .startval(y, Z, tol)
- phi <- startval$phi
- s[p] <- startval$s
- if (equal.LS) { # for testing purpose only
- psi1 <- "Ident"
- niter <- maxiter
- IWLS <- .IWLS(y, Z, phi, s[p], u, v, psi1, niter, tol)
- phi <- IWLS$phi
- w <- IWLS$w
- niter.testing <- IWLS$niter
- } else {
- if ((iterh > 0) & (is.numeric(iterh))) {
- psi1 <- "Huber"
- niter <- iterh
- IWLS <- .IWLS(y, Z, phi, s[p], u, v, psi1, niter, tol, k=chr)
- phi <- IWLS$phi
- s[p] <- IWLS$s
- w <- IWLS$w
- BH <- IWLS$B
- niterh <- IWLS$niter
- }
- if ((iterb > 0) & (is.numeric(iterb))) {
- psi1 <- "Tukey"
- niter <- iterb
- IWLS <- .IWLS(y, Z, phi, s[p], u, v, psi1, niter, tol, c=cbr)
- phi <- IWLS$phi
- s[p] <- IWLS$s
- w <- IWLS$w
- BB <- IWLS$B
- niterb <- IWLS$niter
- }
- }
- if (p > 1) {
- Phi <- cbind(rep(0, (p-1)), Phi)
- Phi <- rbind(phi, Phi)
- } else {
- Phi <- as.matrix(phi)
- }
- }
- Cx <- solve(invCp)
-
- return(list(ar=phi, sinnov=s, Cx=Cx, mu=x.huber$mu, sx=sx, u=u, v=v, w=w,
- BH=BH, BB=BB,
- niterh=niterh, niterb=niterb, niter.testing=niter.testing))
-
-}
Deleted: pkg/robKalman/R/arGMinternal.R
===================================================================
--- pkg/robKalman/R/arGMinternal.R 2025-05-15 14:31:23 UTC (rev 82)
+++ pkg/robKalman/R/arGMinternal.R 2025-05-15 14:34:21 UTC (rev 83)
@@ -1,242 +0,0 @@
-.ARmodel <- function (x, p)
-{
-##################################
-##
-## R-function: .ARmodel - creates design matrix 'Z' and
-## response vector 'y' of an AR(p) model (internal function)
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-21)
-##
-##################################
-
-## x ... univarite time series (vector)
-## p ... order of AR(p) process
-
- n <- length(x)
- y <- x[(p+1):n]
- Z <- embed(x, p)[1:(n-p), ]
- return(list(y=y, Z=as.matrix(Z)))
-}
-
-.invCp <- function (p, s, Phi)
-{
-##################################
-##
-## R-function: .invCp - computes the inverse p x p covariance matrix
-## (internal function)
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-26)
-##
-##################################
-
-## p ... order of AR(p) process
-## s ... vector of 'sx' and innovation scale estimates for
-## AR(p-1) models of order 1 to (p-1)
-## Phi ... (p-1)x(p-1) matrix of AR(p-1) model parameters
-
- if (p > 1) {
- M1 <- matrix(rep(rev(s), p), p)
- M2 <- cbind(rep(0, (p-1)), (-1)*Phi)
- M2 <- rbind(M2, rep(0, p))
- diag(M2) <- 1
- Ap <- (1/M1)*M2
- invCp <- t(Ap)%*%Ap
- } else {
- invCp <- 1/s[1]^2
- }
- return(invCp)
-}
-
-.Weights <- function (p, Z, invCp, type, psi2, c)
-{
-##################################
-##
-## R-function: .Weights - computes weights for
-## Mallows- or Schweppe-type GM-estimates (internal function)
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-26)
-##
-##################################
-
-## p ... order
-## Z ... AR(p) model matrix
-## invCp ... martix from function '.invCp' to compute metric
-## type ... Mallows- or Schweppe-type GM-estimates (character)
-## psi2 ... influence function (character)
-## c ... tuning constant
-
- psi <- .psi(psi2)
- d <- sqrt(diag(Z%*%invCp%*%t(Z))/p)
- if (psi2 == "Huber") {
- v <- psi(d, k=c)/d
- } else if (psi2 == "Tukey") {
- v <- psi(d, c=c)/d
- } else if (psi2 == "Ident") {
- v <- rep(1, nrow(Z))
- } else {
- warning("error in function '.Weights': psi function ", psi2,
- "not defined \n")
- }
- if (type=="Mallows") {
- u <- rep(1, length(v))
- } else if (type=="Schweppe") {
- u <- v
- } else {
- warning("error in function '.Weights': wrong GM-estimates type \n")
- }
- return(list(u=u, v=v))
-}
-
-.startval <- function (y, Z, tol)
-{
-##################################
-##
-## R-function: .startval - computes appropriate starting values
-## (internal function)
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-26)
-##
-##################################
-
-## y ... response vector of AR(p) model
-## Z ... design matrix of AR(p) model
-## tol ... tolerance level
-
- type <- FALSE # hard-coded, TRUE as proposed in [MarT82b], [StoDut87]
- # FALSE as proposed in [MarZ78]
- ar.ls <- lm.fit(Z, y, tol)
- phi <- ar.ls$coefficients
- if (type) {
- s <- sqrt(sum((ar.ls$residuals)^2)/ar.ls$df.residual)
- } else {
- s <- mad(ar.ls$residuals)
- }
- return(list(phi=phi, s=s))
-}
-
-.BH <- function (k=1.345)
-{
-##################################
-##
-## R-function: .BH - computes appropriate constant to obtain a
-## consistent estimate for sigma when using Huber's psi function
-## (internal function)
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-26)
-##
-##################################
-
-## k ... tuning constant (default: k=1.345)
-
- -2*k*dnorm(k) + 2*pnorm(k) + 2*k^2*(1-pnorm(k)) -1
-}
-
-.BB <- function (c=4.685)
-{
-##################################
-##
-## R-function: .BB - computes appropriate constant to obtain a
-## consistent estimate for sigma when using Tukey's psi function
-## (internal function)
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-26)
-##
-##################################
-
-## c ... tuning constant (default: c=4.685)
-
- (2/c^7)*dnorm(c)*(c^6 - 13*c^4 + 105*c^2 - 945) +
- (2/c^8)*pnorm(c)*(c^8 - 12*c^6 + 90*c^4 - 420*c^2 + 945) -
- (c^8 - 12*c^6 + 90*c^4 - 420*c^2 + 945)/c^8
-}
-
-.weights <- function (r, s, u, v, psi1, ...)
-{
-##################################
-##
-## R-function: .weights - computes appropriate weights for reweighting
-## (internal function)
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-26)
-##
-##################################
-
-## r ... residuals
-## s ... innovations scale parameter
-## u ... weights
-## v ... weights
-## psi1 ... influence function (character)
-## ... ... passing tuning constants to influence functions
-
- psi <- .psi(psi1)
- n <- length(r)
- w <- rep(NA, n)
- for (i in 1:n) {
- if (r[i] == 0) {
- if (u[i] != 0) {
- w[i] <- v[i]/u[i]
- } else {
- w[i] <- 1
- }
- } else if (u[i] != 0) {
- dummy <- r[i]/s
- w[i] <- v[i]*psi(dummy/u[i], ...)/dummy
- } else if (psi1 == "Ident") {
- w[i] <- 1
- } else {
- w[i] <- 0
- }
- }
- return(w)
-}
-
-.IWLS <- function (y, Z, phi.ini, s.ini, u, v, psi1, niter, tol, ...)
-{
-##################################
-##
-## R-function: .IWLS - iteratively reweighted least squares algorithm
-## (internal function)
-## author: Bernhard Spangl
-## version: 1.0 (2006-05-26)
-##
-##################################
-
-## y ... response vector of AR(p) model
-## Z ... design matrix of AR(p) model
-## phi.ini ... initial AR(p) model parameters
-## s.ini ... initial innovations scale parameter
-## u ... weights
-## v ... weights
-## psi1 ... influence function (character)
-## niter ... maximal number of iterations
-## tol ... tolerance level
-## ... ... passing tuning constants to influence functions
-
- stop1 <- sqrt(diag(solve(t(Z)%*%Z)))
- B <- NA
- phi <- phi.ini
- s.new <- s.ini
- iter <- 0
- psi <- .psi(psi1)
- if (psi1=="Huber") {
- B <- .BH(...)
- } else {
- B <- .BB(...)
- }
-
- while (iter < niter) {
- iter <- iter + 1
- r <- y - Z%*%phi
- s.old <- s.new
- s2 <- (1/(B*sum(u*v)))*sum(u*v*(psi(r/(u*s.old), ...))^2)*s.old^2
- s.new <- sqrt(s2)
- w <- .weights(r, s.new, u, v, psi1, ...)
- ar.wls <- lm.wfit(Z, y, w, tol)
- tau <- ar.wls$coefficients - phi
- omega <- 1 # hard-coded, 0 < omega < 2, as proposed in [StoDut87]
- phi <- phi + omega*tau
- stop2 <- abs(c(s.old - s.new, omega*tau)) < tol*s.new*c(1, stop1)
- if (sum(stop2) == length(stop2)) break
- }
- return(list(phi=phi, s=s.new, w=w, B=B, niter=iter))
-}
More information about the Robkalman-commits
mailing list