Learn R Programming
Rfssa (version 2.1.0)

fssa: Functional Singular Spectrum Analysis

Description

This is a function which performs the decomposition (including embedding and functional SVD steps) stage for univariate functional singular spectrum analysis (ufssa) or multivariate functional singular spectrum analysis (mfssa). The algorithm (ufssa or mfssa) is chosen based on whether the supplied input is a univariate or multivariate functional time series (fts) object. The type parameter can also be set to "mfssa" if the user wishes to perform ufssa of a univariate fts object using the mfssa code. Also note that the variables of the fts maybe observed over different dimensional domains where the maximum dimension currently supported is two.

Usage

fssa(Y, L = NA, ntriples = 20, type = "ufssa")

Value

An object of class fssa, which is a list of functional objects and the following components:

values

A numeric vector of eigenvalues.

L

The specified window length.

N

The length of the functional time series.

Y

The original functional time series.

Arguments

Y

An object of class fts.

L

A positive integer giving the window length.

ntriples

A positive integer specifying the number of eigentriples to calculate in the decomposition.

type

A string indicating which type of fssa to perform. Use type="ufssa" to perform univariate fssa (default for univariate fts). Use type="mfssa" to perform multivariate fssa (default for multivariate fts).

Examples

Run this code
if (FALSE) {
## Univariate FSSA Example on Callcenter data
require(Rfssa)
load_github_data("https://github.com/haghbinh/Rfssa/blob/master/data/Callcenter.RData")
## Define functional objects
D <- matrix(sqrt(Callcenter$calls), nrow = 240)
N <- ncol(D)
time <- substr(seq(ISOdate(1999, 1, 1), ISOdate(1999, 12, 31), by = "day"),1,10)
K <- nrow(D)
u <- seq(0, K, length.out = K)
d <- 22 # Optimal Number of basis elements
## Define functional time series
Y <- Rfssa::fts(list(D), list(list(d, "bspline")), list(u),time)
Y
plot(Y, mains = c("Call Center Data Line Plot"),
xlabels = "Time (6 minutes aggregated)",
ylabels = "Sqrt of Call Numbers",type="line",
xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
 list(c(1,60,120,180,240)))
## Univariate functional singular spectrum analysis
L <- 28
U <- fssa(Y, L)
plot(U, d = 13)
plot(U, d = 9, type = "lheats")
plot(U, d = 9, type = "lcurves")
plot(U, d = 9, type = "vectors")
plot(U, d = 10, type = "periodogram")
plot(U, d = 10, type = "paired")
plot(U, d = 10, type = "wcor")
gr <- list(1, 2:3, 4:5, 6:7, 1:7)
Q <- freconstruct(U, gr)
plot(Q[[1]], mains = "Call Center Mean Component",
xlabels = "Time (6 minutes aggregated)",
     ylabels = "Sqrt of Call Numbers",type="line",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))
plot(Q[[2]], mains = "Call Center First Periodic Component",
xlabels = "Time (6 minutes aggregated)",
     ylabels = "Sqrt of Call Numbers",type="line",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))
plot(Q[[3]], mains = "Call Center Second Periodic Component",
xlabels = "Time (6 minutes aggregated)",
     ylabels = "Sqrt of Call Numbers",type="line",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))
plot(Q[[4]], mains = "Call Center Third Periodic Component",
xlabels = "Time (6 minutes aggregated)",
     ylabels = "Sqrt of Call Numbers",type="line",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))
plot(Y, mains = c("Call Center Data Line Plot"),
xlabels = "Time (6 minutes aggregated)",
     ylabels = "Sqrt of Call Numbers",type="line",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))
plot(Q[[5]], mains = "Call Center Extracted Signal",
xlabels = "Time (6 minutes aggregated)",
     ylabels = "Sqrt of Call Numbers",type="line",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))

## Other visualization types for object of class "fts":

plot(Q[[1]],xlabels = "Time (6 minutes aggregated)",
     zlabels = "Sqrt of Call Numbers",type="3Dsurface", tlabels = "Date",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))
plot(Q[[2]], mains = "Call Center First Periodic Component",
xlabels = "Time (6 minutes aggregated)",
     zlabels = "Sqrt of Call Numbers",type="heatmap",
     tlabels = "Date",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))
plot(Q[[3]],xlabels = "Time (6 minutes aggregated)",
     zlabels = "Sqrt of Call Numbers",type="3Dline",
     tlabels = "Date",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00")),xticklocs =
       list(c(1,60,120,180,240)))

## Multivariate FSSA Example on bivariate intraday
## temperature curves and smoothed images of vegetation
require(Rfssa)
load_github_data("https://github.com/haghbinh/Rfssa/blob/master/data/Montana.RData")
Temp <- Montana$Temp
NDVI <- Montana$NDVI
d_temp <- 11
d_NDVI <- 13
## Define functional time series
Y <- Rfssa::fts(
  list(Temp / sd(Temp), NDVI), list(
    list(d_temp, "bspline"),
    list(d_NDVI, d_NDVI, "bspline", "bspline")
  ),
  list(c(0, 23), list(c(1, 33), c(1, 33))),
colnames(Temp))
# Plot the first 100 observations
plot(Y[1:100],
     xlabels = c("Time", "Longitude"),
     ylabels = c("Normalized Temperature (\u00B0C)", "Latitude"),
     zlabels = c("", "NDVI"),
     mains = c("Temperature Curves", "NDVI Images"), color_palette = "RdYlGn",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00"),
     c("113.40\u00B0 W", "113.30\u00B0 W")),xticklocs =
       list(c(1,6,12,18,24),c(1,33)),
       yticklabels = list(NA,c("48.70\u00B0 N", "48.77\u00B0 N")),yticklocs =
       list(NA,c(1,33))
)
plot(Y, types = c("3Dline", "heatmap"), vars = c(1, 1))
plot(Y, types = "heatmap", vars = 2)
plot(Y, vars = c(2, 1))
L <- 45
## Multivariate functional singular spectrum analysis
U <- fssa(Y, L)
plot(U, type = "values", d = 10)
plot(U, type = "vectors", d = 4)
plot(U, type = "lheats", d = 4)
plot(U, type = "lcurves", d = 4, vars = c(1))
plot(U, type = "paired", d = 6)
plot(U, type = "wcor", d = 10)
plot(U, type = "periodogram", d = 4)
# Reconstruction of multivariate fts observed over different dimensional domains
Q <- freconstruct(U = U, groups = list(c(1), c(2:3), c(4)))
# Plotting reconstructions to show accuracy
plot(Q[[1]],
     xlabels = c("Time", "Longitude"),
     ylabels = c("Normalized Temperature (\u00B0C)", "Latitude"),
     zlabels = c("", "NDVI"),
     mains = c("Temperature Curves Mean", "NDVI Images Mean"), color_palette = "RdYlGn",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00"),
     c("113.40\u00B0 W", "113.30\u00B0 W")),xticklocs =
       list(c(1,6,12,18,24),c(1,33)),
       yticklabels = list(NA,c("48.70\u00B0 N", "48.77\u00B0 N")),yticklocs =
       list(NA,c(1,33))) # mean
plot(Q[[2]],
     xlabels = c("Time", "Longitude"),
     ylabels = c("Normalized Temperature (\u00B0C)", "Latitude"),
     zlabels = c("", "NDVI"),
     mains = c("Temperature Curves Periodic", "NDVI Images Periodic"), color_palette = "RdYlGn",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00"),
     c("113.40\u00B0 W", "113.30\u00B0 W")),xticklocs =
       list(c(1,6,12,18,24),c(1,33)),
       yticklabels = list(NA,c("48.70\u00B0 N", "48.77\u00B0 N")),yticklocs =
       list(NA,c(1,33))) # periodic
plot(Q[[3]],
     xlabels = c("Time", "Longitude"),
     ylabels = c("Normalized Temperature (\u00B0C)", "Latitude"),
     zlabels = c("", "NDVI"),
     mains = c("Temperature Curves Trend", "NDVI Images Trend"), color_palette = "RdYlGn",
     xticklabels = list(c("00:00","06:00","12:00","18:00","24:00"),
     c("113.40\u00B0 W", "113.30\u00B0 W")),xticklocs =
       list(c(1,6,12,18,24),c(1,33)),
       yticklabels = list(NA,c("48.70\u00B0 N", "48.77\u00B0 N")),yticklocs =
       list(NA,c(1,33))) # trend
}

Run the code above in your browser using DataLab