[Robust-ts-commits] r20 - in pkg/robust-ts: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 23 18:51:54 CET 2015
Author: ruckdeschel
Date: 2015-02-23 18:51:54 +0100 (Mon, 23 Feb 2015)
New Revision: 20
Added:
pkg/robust-ts/R/arGM.R
pkg/robust-ts/R/arGMinternal.R
pkg/robust-ts/man/arGM.Rd
pkg/robust-ts/man/internalarGM.Rd
Log:
a long forgotten commit... Bernhard Spangl's arGM routines
Added: pkg/robust-ts/R/arGM.R
===================================================================
--- pkg/robust-ts/R/arGM.R (rev 0)
+++ pkg/robust-ts/R/arGM.R 2015-02-23 17:51:54 UTC (rev 20)
@@ -0,0 +1,114 @@
+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))
+
+}
Added: pkg/robust-ts/R/arGMinternal.R
===================================================================
--- pkg/robust-ts/R/arGMinternal.R (rev 0)
+++ pkg/robust-ts/R/arGMinternal.R 2015-02-23 17:51:54 UTC (rev 20)
@@ -0,0 +1,242 @@
+.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))
+}
Added: pkg/robust-ts/man/arGM.Rd
===================================================================
--- pkg/robust-ts/man/arGM.Rd (rev 0)
+++ pkg/robust-ts/man/arGM.Rd 2015-02-23 17:51:54 UTC (rev 20)
@@ -0,0 +1,71 @@
+\name{arGM}
+\alias{arGM}
+
+\title{GM-estimates for AR parameters}
+
+\description{
+This function realizes Martin's GM-estimates for AR parameters.
+}
+
+\usage{
+arGM(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, ...)
+}
+
+\arguments{
+\item{x}{univarite time series (vector)}
+\item{order}{order of AR(p) process}
+\item{chr}{tuning constant for Huber's \eqn{\psi} function}
+\item{iterh}{number of iterations for IWLS-alogrithm using Huber's \eqn{\psi} function}
+\item{cbr}{tuning constant for Tukey's \eqn{\psi} function}
+\item{iterb}{number of iterations for IWLS-alogrithm using Tukey's \eqn{\psi} function}
+\item{psi2}{influence function to determine the 'size of \eqn{z_i}', either "Ident", "Huber" or "Tukey"}
+\item{c}{tuning constant for psi2}
+\item{type}{type of GM-estimates, either "Mallows" or "Schweppe"}
+\item{k}{tuning constant for centering}
+\item{maxiter}{maximal number of iteration}
+\item{tol}{tolerance level}
+\item{equal.LS}{logical, for testing purpose only}
+\item{...}{further parameters to be passed to the functions \code{HuberM} or \code{hubers}}
+}
+
+\details{
+to be filled}
+
+\value{a list with elements
+\item{ar}{parameter estimate}
+\item{sinnov}{scale estimate for the innovations from AR(p) fits of orders 1 through p}
+\item{Cx}{an estimate of the \eqn{p \times p}{p x p} autocovariance matrix}
+\item{mu}{ location estimate of \code{x}}
+\item{sx}{Huber scale estimate}
+\item{u,v}{weights for Mallows- or Schweppe-type GM-estimates}
+\item{w}{weights from IWLS algorithm}
+\item{BH}{consistency constant for \eqn{\sigma} when using Huber's \eqn{\psi} function}
+\item{BB}{consistency constant for \eqn{\sigma} when using Tukey's \eqn{\psi} function}
+\item{niterh}{number of iterations for IWLS-alogrithm using Huber's \eqn{\psi} function}
+\item{niterb}{number of iterations for IWLS-alogrithm using Tukey's \eqn{\psi} function}
+\item{niter.testing}{for testing purposes only}
+}
+
+\author{
+Bernhard Spangl \email{bernhard.spangl at boku.ac.at},\cr
+ }
+
+
+\seealso{
+\code{\link{internalarGM}}, \code{\link{internalpsi}}
+}
+
+\examples{
+}
+
+\references{
+Martin, R.D. and Zeh, J.E. (1978): Generalized M-estimates for Autoregression Including Small-sample Efficiency Robustness \cr
+Martin, R.D. (1980): Robust Estimation of Autoregressive Models. \cr
+Martin, R.D. and Thomson, D.J. (1982): Robust-resistent Spectrum Estimation. \cr
+Stockinger, N. and Dutter, R. (1987): Robust Time Series Analysis: A Survey. \cr
+}
+
+\keyword{robust}
+\keyword{ts}
Added: pkg/robust-ts/man/internalarGM.Rd
===================================================================
--- pkg/robust-ts/man/internalarGM.Rd (rev 0)
+++ pkg/robust-ts/man/internalarGM.Rd 2015-02-23 17:51:54 UTC (rev 20)
@@ -0,0 +1,77 @@
+\name{internalarGM}
+\alias{internalarGM}
+\alias{.ARmodel}
+\alias{.invCp}
+\alias{.Weights}
+\alias{.startval}
+\alias{.BH}
+\alias{.BB}
+\alias{.weights}
+\alias{.IWLS}
+
+\title{Internal functions of package robKalman --- psi functions}
+
+\description{
+These functions are used internally by package robKalman for the ACM filter --- \eqn{\psi}-functions
+}
+
+\usage{
+.ARmodel(x, p)
+.invCp(p, s, Phi)
+.Weights(p, Z, invCp, type, psi2, c)
+.startval(y, Z, tol)
+.BH(k=1.345)
+.BB(c=4.685)
+.weights(r, s, u, v, psi1, ...)
+.IWLS(y, Z, phi.ini, s.ini, u, v, psi1, niter, tol, ...)
+}
+
+\arguments{
+ \item{x}{univarite time series (vector)}
+ \item{p}{order of AR(p) process}
+ \item{s}{vector of \code{sx} and innovation scale estimates for AR(p-1) models of order \eqn{1} to \eqn{(p-1)}}
+ \item{Phi}{\eqn{(p-1)\times(p-1)}{(p-1)x(p-1)} matrix of AR(p-1) model parameters}
+ \item{Z}{AR(p) model matrix}
+ \item{invCp}{matrix from function \code{.invCp} to compute metric}
+ \item{type}{type of GM-estimates --- character: currently: "Mallows" or "Schweppe"}
+ \item{psi1, psi2}{type of \eqn{\psi} function; current possibilities: "Huber", "Tukey", "Hampel", "Ident"}
+ \item{c, k}{tuning constants}
+ \item{y}{response vector of AR(p) model}
+ \item{tol}{tolerance level}
+ \item{r}{residuals}
+ \item{s}{innovations scale parameter}
+ \item{u, v}{weights}
+ \item{...}{additional arguments (tuning constants) for influence functions}
+ \item{phi.ini}{initial AR(p) model parameters}
+ \item{s.ini}{initial innovations scale parameter}
+ \item{niter}{maximal number of iterations}
+}
+
+\details{
+to be filled
+}
+
+
+\value{
+\code{.ARmodel(x, p)} returns design matrix \code{Z} and response vector \code{y} of an AR(p) model (as list with corresponding elements).\cr
+\code{.invCp(p, s, Phi)} computes the inverse \eqn{p \times p}{p x p} covariance matrix.\cr
+\code{.Weights} computes weights for Mallows- or Schweppe-type GM-estimates; returns a list with elements \code{u}, \code{v}.\cr
+\code{.startval} computes appropriate starting values; returns a list with elements \code{phi}, \code{s}.\cr
+\code{.BH} computes appropriate constant to obtain a consistent estimate for \eqn{\sigma} when using Huber's \eqn{\psi} function.\cr
+\code{.BB} computes appropriate constant to obtain a consistent estimate for \eqn{\sigma} when using Tukey's \eqn{\psi} function.\cr
+\code{.weights} computes appropriate weights for reweighting --- returns a vector.\cr
+\code{.IWLS} iteratively reweighted least squares algorithm; returns a list with elements
+ \code{phi}, \code{s}, \code{w}, \code{B}, \code{niter}.\cr
+}
+
+\author{
+Bernhard Spangl \email{bernhard.spangl at boku.ac.at},\cr
+}
+
+
+\seealso{
+\code{\link{internalpsi}}
+}
+
+
+\keyword{internal}
More information about the Robust-ts-commits
mailing list