[Robust-ts-commits] r9 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 4 20:25:25 CET 2009
Author: bspangl
Date: 2009-03-04 20:25:25 +0100 (Wed, 04 Mar 2009)
New Revision: 9
Added:
pkg/R/ACF.R
pkg/R/CovQn.R
pkg/R/PACF.R
pkg/R/PSD3.R
pkg/R/Spearman.R
Log:
Upload of R-functions concerning robust autocovariance estimation, cf. Spangl, 2008, PhD-thesis, Chp. 4 - Bernhard Spangl, 2009-03-04
Added: pkg/R/ACF.R
===================================================================
--- pkg/R/ACF.R (rev 0)
+++ pkg/R/ACF.R 2009-03-04 19:25:25 UTC (rev 9)
@@ -0,0 +1,59 @@
+ACF <- function (x, type="spearman", cov=FALSE, modified=TRUE)
+{
+######################################################
+##
+## R-function 'ACF'
+## computes the autocorrelation function/sequence
+## using Spearman's rank correlation coefficient
+## or highly robust autocovariance function proposed
+## by Ma & Genton (2000),
+## (non-parametric alternative to the autocorrelation
+## coefficient)
+##
+## author: Bernhard Spangl
+## date: 2008-06-01, 2007-04-06
+##
+## x ... time series (vector)
+## type ... which correlation coefficient should be used,
+## Spearman's rho ("spearman") or Qn proposed by
+## Rousseeuw & Croux (1993) ("Qn"),
+## default = "spearman"
+## cov ... should the autocovariance function also
+## be computed in case of using 'type="Qn"'
+## modified ... should the modified autocovariance function
+## be computed in case of using 'type="Qn"'
+##
+######################################################
+
+ n <- length(x)
+ mlag <- n-2
+ acovf <- rep(NA, mlag+1)
+ acf <- rep(NA, mlag+1)
+
+## cor <- switch(type, spearman=get("Spearman", mode="function"),
+## Qn=get("CovQn", mode="function"))
+
+ for (i in 0:mlag) {
+ x1 <- x[(1+i):n]
+ x2 <- x[1:(n-i)]
+ if (type == "spearman") { # Spearman's rho
+ acf[i+1] <- Spearman(x1, x2, n)
+ }
+ if (type == "Qn") { # highly robust scale estimator 'Qn'
+ if (modified) {
+ res <- CovQn(x1, x2, cov, n)
+ acf[i+1] <- res$acf.m
+ if (cov) {
+ acovf[i+1] <- res$acovf.m
+ }
+ } else {
+ res <- CovQn(x1, x2, cov)
+ acf[i+1] <- res$acf
+ if (cov) {
+ acovf[i+1] <- res$acovf
+ }
+ }
+ }
+ }
+ return(list(ACF=c(acf, 0), ACovF=c(acovf, 0)))
+}
Added: pkg/R/CovQn.R
===================================================================
--- pkg/R/CovQn.R (rev 0)
+++ pkg/R/CovQn.R 2009-03-04 19:25:25 UTC (rev 9)
@@ -0,0 +1,44 @@
+CovQn <- function (u, v, cov=FALSE, N)
+{
+######################################################
+##
+## R-function 'CovQn'
+## computes the covariance and correlation coefficient
+## using the highly robust scale estimator 'Qn', to be
+## applied to time series (as non-parametric alternative
+## to the autocovariance and autocorrelation coefficient)
+## (cf. Ma and Genton, 2000)
+##
+## IMPORTANT: requires R-package 'robustbase'
+##
+## author: Bernhard Spangl
+## date: 2008-06-01
+##
+## u, v ... appropriate vectors
+## cov ... boolean, should the autocovariance also be
+## calculated?
+## N ... length of time series for modified estimates
+##
+######################################################
+
+ acf <- acovf <- NA
+ acf.m <- acovf.m <- NA
+ n <- length(u)
+ if (missing(N)) {
+ N <- n
+ }
+ mod.fact <- n/N
+
+ Q1 <- Qn(u+v)^2
+ Q2 <- Qn(u-v)^2
+
+ acf <- (Q1 - Q2)/(Q1 + Q2)
+ acf.m <- mod.fact*acf
+ if (cov) {
+ acovf <- (Q1 - Q2)/4
+ acovf.m <- mod.fact*acovf
+ }
+
+ return(list(acf=acf, acf.m=acf.m, acovf=acovf, acovf.m=acovf.m))
+
+}
Added: pkg/R/PACF.R
===================================================================
--- pkg/R/PACF.R (rev 0)
+++ pkg/R/PACF.R 2009-03-04 19:25:25 UTC (rev 9)
@@ -0,0 +1,57 @@
+AcfPacf <- function (x, mx, type="spearman")
+{
+######################################################
+##
+## R-function 'AcfPacf':
+## computes the autocorrelation cofficients recursively
+## via the partial autocorrelation cofficients
+## (cf. Möttönen et al., 1999)
+##
+## author: Bernhard Spangl
+## date: 2007-04-09
+##
+## x ... time series (vector)
+## mx ... mean of time series x
+## type ... which correlation coefficient should be used,
+## Pearson's or Spearman's? (default = spearman)
+##
+######################################################
+
+ if (missing(mx)) {
+ mx <- mean(x)
+ }
+ x <- x - mx # centering
+ n <- length(x)
+ mlag <- n - 2 # maximum lag using this recursive algorithm
+ rho <- rep(NA, n)
+ A <- matrix(NA, n-2, n-2)
+
+ ## Doesn't work for 'pearson', because the biased estimator for
+ ## the Pearson's correlation coefficient has to be used (i.e.,
+ ## should always be divided by length of series)!
+ cor <- switch(type, pearson=get("cor", mode="function"),
+ spearman=get("Spearman", mode="function"))
+
+ rho[1] <- 1 # lag zero
+
+ A[1, 1] <- cor(x[1:(n-1)], x[2:n], n)
+## A[1, 1] <- pacf(x, lag.max=n-1, plot=F)$acf[1, 1, 1] # built-in
+ rho[2] <- A[1, 1]
+
+ for (h in 2:mlag) {
+ X <- embed(x, h+1)
+ u <- X[, 1] - X[, 2:h]%*%t(A[h-1, 1:(h-1), drop=F])
+ v <- X[, h+1] - X[, h:2]%*%t(A[h-1, 1:(h-1), drop=F])
+ A[h, h] <- cor(u, v, n)
+## A[h, h] <- pacf(x, lag.max=n-1, plot=F)$acf[h, 1, 1] # built-in
+ for (i in 1:(h-1)) {
+ A[h, i] <- A[h-1, i] - A[h, h]*A[h-1, h-i]
+ }
+ rho[h+1] <- A[h-1, 1:(h-1)]%*%rev(rho[2:h]) +
+ A[h, h]*(1 - A[h-1, 1:(h-1)]%*%rho[2:h])
+ }
+
+ rho[n] <- 0
+ return(list(ACF=rho, A=A))
+
+}
Added: pkg/R/PSD3.R
===================================================================
--- pkg/R/PSD3.R (rev 0)
+++ pkg/R/PSD3.R 2009-03-04 19:25:25 UTC (rev 9)
@@ -0,0 +1,59 @@
+PSD3 <- function (r, f, finv, delta=0.05, ...)
+{
+######################################################
+##
+## R-function 'PSD3' - non-linear shrinking:
+## computes a positiv semidefinite autocorrelation function
+## via the shrinking method (cf. Rousseeuw and Molenberghs, 1993),
+## using implementations (function 'is.positive.definite')
+## of R package 'corpcor' (cf. Higham, 1988),
+## cf. also Devlin et al. (1975)
+##
+## author: Bernhard Spangl
+## date: 2008-06-11
+##
+## r ... originial autocorrelation function
+## f ... monotone increasing continuous function
+## finv ... inverse of 'f'
+## delta ... shifting parameter (e.g.: delta=0.05)
+## ... ... tolerance (e.g.: tol=1e-08)
+##
+######################################################
+
+ r <- as.vector(r)
+ R <- toeplitz(r)
+
+ was.psd <- is.positive.definite(R, method="eigen", ...)
+ is.psd <- was.psd
+
+ if (!was.psd) {
+ limit <- finv(delta)
+ while (!is.psd) {
+ r <- sapply(r[-1], Trans, f=f, finv=finv, delta=delta, limit=limit)
+ r <- c(1, r)
+ R.psd <- toeplitz(r)
+ is.psd <- is.positive.definite(R.psd, method="eigen", ...)
+ }
+ res <- list(r.psd=R.psd[1, ], R.psd=R.psd, R=R,
+ was.psd=was.psd, is.psd=is.psd)
+ } else {
+ res <- list(r.psd=r, R.psd=R, R=R,
+ was.psd=was.psd, is.psd=is.psd)
+ }
+
+ return(res)
+
+}
+
+Trans <- function (r, f, finv, delta, limit)
+{
+ if (r < -limit) {
+ finv(f(r) + delta)
+ } else {
+ if (r <= limit) {
+ 0
+ } else {
+ finv(f(r) - delta)
+ }
+ }
+}
Added: pkg/R/Spearman.R
===================================================================
--- pkg/R/Spearman.R (rev 0)
+++ pkg/R/Spearman.R 2009-03-04 19:25:25 UTC (rev 9)
@@ -0,0 +1,29 @@
+Spearman <- function (x, y, N)
+{
+######################################################
+##
+## R-function 'Spearman'
+## computes Spearman's (modified) rank correlation
+## coefficient, to be applied to time series (as
+## non-parametric alternative to the autocorrelation
+## coefficient)
+## (cf. Ahdesmäki et al., 2005, Equation 12)
+##
+## author: Bernhard Spangl
+## date: 2007-04-06
+##
+## x, y ... vectors
+## N ... interger, should biased or unbiased version be used?
+##
+######################################################
+
+ if (missing(N)) {
+ N <- length(x)
+ }
+ n <- length(x)
+ rx <- rank(x)
+ ry <- rank(y)
+## rho <- cor(x, y, method="spear")*n/N # build-in version
+ rho <- 12/(N*(n^2 - 1))*((rx - (n + 1)/2)%*%(ry - (n + 1)/2))
+ return(drop(rho))
+}
More information about the Robust-ts-commits
mailing list