Nothing
misst <- function(y, l = 'automatic', m = 'automatic', vertical = TRUE)
UseMethod("misst")
misst.default <- function(y, l = 'automatic', m = 'automatic',
vertical = TRUE) {
n <- y$n
D <- y$D
## H. Hassani & R. Mahmoudvand (2013; pp. 68, eq. 25).
if((l == 'automatic') & (vertical == T | vertical == TRUE)) {l <- ceiling((n + 1) / (D + 1))}
if((l == 'automatic') & (vertical == F | vertical == FALSE)){l <- ceiling(D * (n + 1) / (D + 1))}
k <- n - l + 1
## Run a basic input validation
if (class(y) != "mitsframe")
stop('Stop: the input y must be a multivariate interval valued object (see ?mitsframe).')
if (is.numeric(l) & (l%%1 != 0 | l <= 0 | l > y$n))
stop('Stop: l must be a positive integer smaller than number of observations per series.')
if(is.numeric(m) & length(m) != D)
stop('Stop: The length of vector m must be equal to the number of time series to analyze.')
## Step 1:
Y_list = list();
for(d in 1:D){
Y <- array(NA, c(l, k, 2) ) # l rows, k columns, 2 dimensions
Y[ , , 1] <- trajectory(y$A[, d], l, k)
Y[ , , 2] <- trajectory(y$B[, d], l, k)
Y_list[[d]] <- Y
}
if(vertical == T | vertical == TRUE) {
Y <- array(NA, c(D*l, k,2) )
for( i in 1:D){ Y[((i-1)*l+1):(i*l) , ,] <- Y_list[[i]] } } else {
Y <- array(NA, c(l, D*k,2) )
for( i in 1:D){ Y[ ,((i-1)*k+1):(i*k), ] <- Y_list[[i]] } }
## Step 2: SVD
if(vertical == T | vertical == TRUE) {
S <- matrix(NA, nrow = D*l, ncol = D*l)
for( i in 1:D){
for(j in 1:D){
S[((i-1)*l+1):(i*l),((j-1)*l+1):(j*l)] <- Var.Est(Y_list[[i]],Y_list[[j]])$outer.product;
}
} } else {
S <- matrix(0, nrow = l, ncol = l)
for( i in 1:D){
S = S + Var.Est(Y_list[[i]],Y_list[[i]])$outer.product
}
} # isSymmetric(S)
SVD <- svd(S)
rk <- qr(SVD$u)$rank # rank.
## Step 3 and 4:
A.tilde = B.tilde = A.Residuals = B.Residuals = matrix(0, ncol = D, nrow = n)
if(is.numeric(m) == T) {
if(sum(m > rk) > 0) {
m = replace(m, m > rk, rk)
warning( "Number of components (see 'm') required to construct trendlines in at least
^ one dimension was reduced to be compatible with the rank of the trajectory matrix." ) }
if(vertical == TRUE){
for(d in 1:D){
Y_i = YYI_i(Y, SVD, i=m[d], l=D*l, k=k) ;
tempA <- dbar(Y_i$Y_iA[((d - 1) * l + 1): (d * l), ], l = l, k = k)
tempB <- dbar(Y_i$Y_iB[((d - 1) * l + 1): (d * l), ], l = l, k = k)
A.tilde[,d] <- pmin(tempA, tempB)
B.tilde[,d] <- pmax(tempA, tempB)
A.Residuals[,d] <- pmin(y$A[,d]-A.tilde[,d],y$B[,d]-B.tilde[,d])
B.Residuals[,d] <- pmax(y$A[,d]-A.tilde[,d],y$B[,d]-B.tilde[,d])
} } else {
for(d in 1:D){
Y_i = YYI_i(Y, SVD, i=m[d], l=l, k=D*k) ;
tempA <- dbar(Y_i$Y_iA[ , ((d - 1) * k + 1) : (d * k) ], l = l, k = k)
tempB <- dbar(Y_i$Y_iB[ , ((d - 1) * k + 1) : (d * k) ], l = l, k = k)
A.tilde[,d] <- pmin(tempA, tempB)
B.tilde[,d] <- pmax(tempA, tempB)
A.Residuals[,d] <- pmin(y$A[,d]-A.tilde[,d],y$B[,d]-B.tilde[,d])
B.Residuals[,d] <- pmax(y$A[,d]-A.tilde[,d],y$B[,d]-B.tilde[,d])
} } } else{
m = rep(0,D); #* Storing here the number of ERC on each dimension.
for(d in 1:D){
stop.flag = 1
while(stop.flag != 0 & m[d] < rk) {
m[d] = m[d] + 1
if(vertical == TRUE){
Y_i = YYI_i(Y, SVD, i=m[d], l=D*l, k=k) ;
tempA <- dbar(Y_i$Y_iA[((d - 1) * l + 1): (d * l), ], l = l, k = k)
tempB <- dbar(Y_i$Y_iB[((d - 1) * l + 1): (d * l), ], l = l, k = k)
A.tilde[,d] <- pmin(tempA, tempB)
B.tilde[,d] <- pmax(tempA, tempB)
A.Residuals[,d] <- pmin(y$A[,d]-A.tilde[,d],y$B[,d]-B.tilde[,d])
B.Residuals[,d] <- pmax(y$A[,d]-A.tilde[,d],y$B[,d]-B.tilde[,d])
} else {
Y_i = YYI_i(Y, SVD, i=m[d], l=l, k=D*k) ;
tempA <- dbar(Y_i$Y_iA[ , ((d - 1) * k + 1) : (d * k) ], l = l, k = k)
tempB <- dbar(Y_i$Y_iB[ , ((d - 1) * k + 1) : (d * k) ], l = l, k = k)
A.tilde[,d] <- pmin(tempA, tempB)
B.tilde[,d] <- pmax(tempA, tempB)
A.Residuals[,d] <- pmin(y$A[,d]-A.tilde[,d],y$B[,d]-B.tilde[,d])
B.Residuals[,d] <- pmax(y$A[,d]-A.tilde[,d],y$B[,d]-B.tilde[,d])
}
stop.flag <- ecip(as.matrix(cbind(A.Residuals[,d],B.Residuals[,d])));
} } } # end 'else'
#### ERC:
erc = list()
for(j in 1:D){
erc.d <- array(NA, c(n, m[j], 2) ) # n rows, m columns, 2 dimensions
for(i in 1:m[j]){
if(vertical == TRUE){
Y_i = YYI_i(Y, SVD, i=i, l=D*l, k=k) ;
ll <- dbar(Y_i$Y_iA[((d - 1) * l + 1): (d * l), ], l = l, k = k)
uu <- dbar(Y_i$Y_iB[((d - 1) * l + 1): (d * l), ], l = l, k = k)
erc.d[,i,1] <- pmin(ll,uu); erc.d[,i,2] <- pmax(ll,uu)} else {
Y_i = YYI_i(Y, SVD, i=i, l=l, k=D*k) ;
ll <- dbar(Y_i$Y_iA[ , ((d - 1) * k + 1) : (d * k) ], l = l, k = k)
uu <- dbar(Y_i$Y_iB[ , ((d - 1) * k + 1) : (d * k) ], l = l, k = k)
erc.d[,i,1] <- pmin(ll,uu); erc.d[,i,2] <- pmax(ll,uu) } }
erc[[j]] <- erc.d
}
# Grouping the results into mtsframe objects to deliver in the output
Y.tilde <- mitsframe(dates = y$dates, A = A.tilde, B = B.tilde)
residuals <- mitsframe(dates = y$dates, A = A.Residuals, B = B.Residuals)
outputs <- list(trendlines = Y.tilde,
l = l,
m = m,
vertical = vertical,
residuals = residuals,
svd = SVD,
erc = erc,
observations = y,
call = match.call()) ;
class(outputs) <- "misst"
return(outputs)
}
print.misst <- function(x, ...) {
cat("\n Multivariate Singular Spectrum Trendlines for Interval Data:\n ========================================= \n")
print(x$call)
cat("\n Interval Trendlines:\n"); cat("\n")
print(x$trendlines)
cat("\n Elementary components included in the estimation \n"); cat("\n")
print(x$m)
}
plot.misst <- function(x, time.format = "%m-%y", col = NULL, lty = NULL, main = NULL, type = NULL, pch = NULL, lwd = NULL,
ylab = NULL,xlab = NULL, ylim = NULL, xlim = NULL, cex.lab = NULL, cex.axis = NULL, cex.main = NULL,
tick = FALSE, options = list(type = 'trendlines' , ncomp = 1:5), ...) {
if(options$type %!in% c('trendlines', 'screeplots', 'components', 'cpgrams')) {
print('options must be one of the strings: "trendlines", "components", "cpgrams", or "screeplots".')}
if(options$type == 'trendlines') {
matplot(x$trendlines$dates, x$trendlines$A,
main = if(is.null(main)){''} else {main},
col = if(is.null(col)){'black'} else {col},
lty = if(is.null(lty)){1} else {lty},
pch = if(is.null(pch)){1} else {pch},
lwd = if(is.null(lwd)){1} else {lwd},
type = if(is.null(type)){'l'} else {type},
xlab = if(is.null(xlab)){'Time'} else {xlab},
ylab = if(is.null(ylab)){'Interval Singular Spectrum Trendline'} else {ylab},
ylim = if(is.null(ylim)){c(min(x$trendlines$A,x$trendlines$B),max(x$trendlines$A,x$trendlines$B))} else {ylim},
xlim = if(is.null(xlim)){c(min(x$trendlines$dates),max(x$trendlines$dates))} else {xlim},
cex.lab = if(is.null(cex.lab)){1} else {cex.lab},
cex.axis = if(is.null(cex.axis)){1} else {cex.axis},
cex.main = if(is.null(cex.main)){1} else {cex.main},
xaxt = "n")
if(inherits(x$trendlines$dates, "Date")==T){ timelabels<-format(x$trendlines$dates,time.format) ;
axis(1,at=x$trendlines$dates,labels=timelabels, tick = tick, cex.axis = if(is.null(cex.axis)){1} else {cex.axis})} else { axis(1,at=x$trendlines$dates, tick = tick, cex.axis = if(is.null(cex.axis)){1} else {cex.axis}) }
for(i in 1:x$trendlines$D){
lines(x$trendlines$dates, x$trendlines$A[,i],
col = if(is.null(col)){'black'} else {col},
lty = if(is.null(lty)){1} else {lty},
pch = if(is.null(pch)){1} else {pch},
lwd = if(is.null(lwd)){1} else {lwd},
type = if(is.null(type)){'l'} else {type})
lines(x$trendlines$dates, x$trendlines$B[,i],
col = if(is.null(col)){'black'} else {col},
lty = if(is.null(lty)){1} else {lty},
pch = if(is.null(pch)){1} else {pch},
lwd = if(is.null(lwd)){1} else {lwd},
type = if(is.null(type)){'l'} else {type})
polygon(c(rev(x$trendlines$dates), x$trendlines$dates),
c(rev(x$trendlines$A[,i]), x$trendlines$B[,i]),
col = if(is.null(col)){'lightgray'} else {col}, border = NA)}
}
if(options$type == 'screeplots') {
if(is.null(options$ncomp)){options$ncomp = c(1:3)}
D <- x$trendlines$D
for(i in 1:D) {
plot(log(x$svd$d[options$ncomp]),
main = if(is.null(main)){'Scree--plot'} else {main},
col = if(is.null(col)){'black'} else {col},
lty = if(is.null(lty)){1} else {lty},
pch = if(is.null(pch)){1} else {pch},
lwd = if(is.null(lwd)){1} else {lwd},
type = if(is.null(type)){'b'} else {type},
cex.lab = if(is.null(cex.lab)){1} else {cex.lab},
cex.axis = if(is.null(cex.axis)){1} else {cex.axis},
cex.main = if(is.null(cex.main)){1} else {cex.main},
xlab = if(is.null(xlab)){'Index'} else {xlab},
ylab = if(is.null(ylab)){'Eigenvalues (log-scale)'} else {ylab})
}
}
if(options$type == 'components') {
if(is.null(options$ncomp)){options$ncomp = 1 ; print('ncomp is missing, only the first erc is ploted')}
coldim <- options$ncomp; rowdim <- x$trendlines$D
par(mfrow = c(rowdim, length(coldim)))
for(i in 1:rowdim) {
for(j in coldim) {
plot(x$trendline$dates, x$erc[[i]][, j, 1],
main = if(is.null(main)){paste("Serie: ",i," - ERC:", j)} else {main},
col = if(is.null(col)){'black'} else {col},
lty = if(is.null(lty)){1} else {lty},
pch = if(is.null(pch)){1} else {pch},
lwd = if(is.null(lwd)){1} else {lwd},
type = if(is.null(type)){'l'} else {type},
xlab = if(is.null(xlab)){''} else {xlab},
ylab = if(is.null(ylab)){''} else {ylab},
ylim = if(is.null(ylim)){ c(min(c(x$erc[[i]][, j, 1],x$erc[[1]][, j, 2])),max(c(x$erc[[i]][, j, 1],x$erc[[i]][, j, 2]))) } else {ylim},
cex.lab = if(is.null(cex.lab)){0.8} else {cex.lab},
cex.axis = if(is.null(cex.axis)){0.8} else {cex.axis},
cex.main = if(is.null(cex.main)){1} else {cex.main} )
lines(x$trendline$dates, x$erc[[i]][, j, 2],
col = if(is.null(col)){'black'} else {col},
lty = if(is.null(lty)){1} else {lty},
pch = if(is.null(pch)){1} else {pch},
lwd = if(is.null(lwd)){1} else {lwd},
type = if(is.null(type)){'l'} else {type})
polygon(c(rev(x$trendline$dates), x$trendline$dates),
c(rev(x$erc[[i]][, j, 1]), x$erc[[i]][, j, 2]),
border = NA, col = if(is.null(col)){'lightgray'} else {col})
}
}
}
if(options$type == "cpgrams") {
D <- x$trendlines$D
for(i in 1:D){
ecip(residuos = cbind(x$residuals$A[,i],x$residuals$B[,i]), plot.flag = TRUE)
}
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.