gSOBI: Generalized SOBI
In tsBSS: Blind Source Separation and Supervised Dimension Reduction for Time Series

Description Usage Arguments Details Value Author(s) References See Also Examples

The gSOBI (generalized Second Order Blind Identification) method for the blind source separation (BSS) problem. The method is designed for separating multivariate time series with or without stochastic volatility. The method is a combination of SOBI and vSOBI with G(x) = x^2 as a nonlinearity function.

gSOBI(X, ...)

## Default S3 method:
gSOBI(X, k1 = 1:12, k2 = 1:3, b = 0.9, eps = 1e-06, maxiter = 1000, ordered = FALSE,
      acfk = NULL, original = TRUE, alpha = 0.05, ...)
## S3 method for class 'ts'
gSOBI(X, ...)
## S3 method for class 'xts'
gSOBI(X, ...)
## S3 method for class 'zoo'
gSOBI(X, ...)

`X`	A numeric matrix or a multivariate time series object of class `ts`, `xts` or `zoo`. Missing values are not allowed.
`k1`	A vector of lags for SOBI part. It can be any non-zero positive integer, or a vector consisting of them. Default is `1:12`.
`k2`	A vector of lags for vSOBI part. It can be any non-zero positive integer, or a vector consisting of them. Default is `1:3`.
`b`	The weight for the SOBI part, 1-b for the vSOBI part. Default is 0.9.
`eps`	Convergence tolerance.
`maxiter`	The maximum number of iterations.
`ordered`	Whether to order components according to their volatility. Default is `FALSE`.
`acfk`	A vector of lags to be used in testing the presence of serial autocorrelation. Applicable only if `ordered = TRUE`.
`original`	Whether to return the original components or their residuals based on ARMA fit. Default is `TRUE`, i.e. the original components are returned. Applicable only if `ordered = TRUE`.
`alpha`	Alpha level for linear correlation detection. Default is 0.05.
`...`	Further arguments to be passed to or from methods.

Assume that a p-variate Y with T observations is whitened, i.e. Y = S^(-1/2)*(X_t - (1/T)*sum_t(X_t)), for t = 1, …, T, where S is the sample covariance matrix of X. The algorithm finds an orthogonal matrix U by maximizing

D(U) = b*sum_k_1(D_k_1(U)) + (1 - b)*sum_k_2(D_k_2(U)) ,

where SOBI part

D_k_1 = sum_i(((1/(T - k_1))*sum_t[(u_i' Y_t)*(u_i' Y_(t + k_1))])^2),

and vSOBI part

D_k_2 = sum_i(((1/(T - k_2))*sum_t[(u_i' Y_t)^2*(u_i' Y_(t + k_2))^2] - (1/(T - k_2)^2*sum_t[(u_i' Y_t)^2]*sum_t[(u_i' Y_(t + k_2))^2])^2))

where b is a value between 0 and 1, and i = 1, …, p, k_1 = 1, …, K_1, k_2 = 1, …, K_2 and t = 1, …, T

The algorithm works iteratively starting with diag(p) as an initial value for an orthogonal matrix U = (u_1, u_2, …, u_p)'.

Matrix T_ikj is a partial derivative of D_kj(U), where j = 1, 2, with respect to u_i. Then T_kj = (T_1kj, …, T_pkj)', where p is the number of columns in Y, and T_j = sum(T_kj), for j = 1, 2. Finally T = b*T_1 + (1 - b)*T_2.

The update for the orthogonal matrix U.new = (TT')^(-1/2)*T is calculated at each iteration step. The algorithm stops when

||U.new - U.old||

is less than eps. The final unmixing matrix is then W = U S^(-1/2).

For ordered = TRUE the function orders the sources according to their volatility. First a possible linear autocorrelation is removed using auto.arima. Then a squared autocorrelation test is performed for the sources (or for their residuals, when linear correlation is present). The sources are then put in a decreasing order according to the value of the test statistic of the squared autocorrelation test. For more information, see lbtest.

A list of class 'bssvol', inheriting from class 'bss', containing the following components:

`W`	The estimated unmixing matrix. If `ordered = TRUE`, the rows are ordered according to the order of the components.
`k1`	The vector of the used lags for the SOBI part.
`k2`	The vector of the used lags for the vSOBI part.
`S`	The estimated sources as time series object standardized to have mean 0 and unit variances. If `ordered = TRUE`, then components are ordered according to their volatility. If `original = FALSE`, the sources with linear autocorrelation are replaced by their ARMA residuals.
`MU`	The mean vector of `X`.

If ordered = TRUE, then also the following components included in the list:

`Sraw`	The ordered original estimated sources as time series object standardized to have mean 0 and unit variances. Returned only if `original = FALSE`.
`fits`	The ARMA fits for the components with linear autocorrelation.
`armaeff`	A logical vector. Is TRUE if ARMA fit was done to the corresponding component.
`linTS`	The value of the modified Ljung-Box test statistic for each component.
`linP`	p-value based on the modified Ljung-Box test statistic for each component.
`volTS`	The value of the volatility clustering test statistic.
`volP`	p-value based on the volatility clustering test statistic.

Markus Matilainen, Jari Miettinen

Belouchrani, A., Abed-Meriam, K., Cardoso, J.F. and Moulines, R. (1997), A Blind Source Separation Technique Using Second-Order Statistics, IEEE Transactions on Signal Processing, 434–444.

Matilainen, M., Miettinen, J., Nordhausen, K., Oja, H. and Taskinen, S. (2017), On Independent Component Analysis with Stochastic Volatility Models, Austrian Journal of Statistics, 46(3–4), 57–66.

Miettinen, M., Matilainen, M., Nordhausen, K. and Taskinen, S. (2020), Extracting Conditionally Heteroskedastic Components Using Independent Component Analysis, Journal of Time Series Analysis, 41, 293–311.

SOBI, vSOBI, lbtest, auto.arima

if(require("stochvol")) {
n <- 10000
A <- matrix(rnorm(9), 3, 3)

# simulate SV models
s1 <- svsim(n, mu = -10, phi = 0.8, sigma = 0.1)$y
s2 <- svsim(n, mu = -10, phi = 0.9, sigma = 0.2)$y
s3 <- svsim(n, mu = -10, phi = 0.95, sigma = 0.4)$y

# create a daily time series
X <- ts(cbind(s1, s2, s3) %*% t(A), end = c(2015, 338), frequency = 365.25)

res <- gSOBI(X, 1:4, 1:2, 0.99)
res$W
coef(res)
plot(res)
head(bss.components(res))

MD(res$W, A) # Minimum Distance Index, should be close to zero

# xts series as input
library("xts")
data(sample_matrix)
X2 <- as.xts(sample_matrix)
res2 <- gSOBI(X2, 1:4, 1:2, 0.99)
plot(res2, multi.panel = TRUE)

# zoo series as input
X3 <- as.zoo(X)
res3 <- gSOBI(X3, 1:4, 1:2, 0.99)
plot(res3)
}