swh:1:snp:3a44eb759780145deea094ac2a25c5049546a085
Tip revision: b670a15ed46f9ec07083d3c6c338a20fd3b3e3b9 authored by Han Lin Shang on 21 April 2019, 14:10:06 UTC
version 5.5
version 5.5
Tip revision: b670a15
ftsm.R
ftsm <- function (y, order = 6, ngrid = max(500, ncol(y$y)), method = c("classical",
"M", "rapca"), mean = TRUE, level = FALSE, lambda = 3, weight = FALSE,
beta = 0.1, ...)
{
if(length(colnames(y$y)) > 0)
{
y$time = ts(as.numeric(colnames(y$y)), start = head(as.numeric(colnames(y$y)),1),
end = tail(as.numeric(colnames(y$y)),1), frequency = 1/diff(as.numeric(colnames(y$y)))[1])
}
else
{
y$time = ts(1:ncol(y$y), start=1, end = ncol(y$y))
}
method <- match.arg(method)
if (!mean & !level & order < 1)
stop("No model to fit")
if (order < 0)
stop("Order must not be negative")
n <- ncol(y$y)
m <- length(y$x)
if (weight == FALSE) {
# remove years with less than two data points
yy <-y$y
yynotNA <- nrow(yy) - colSums(is.na(yy)) >= 2
yy <- yy[, yynotNA]
y.pca <- fdpca(y$x, yy, order = order, ngrid = ngrid,
method = method, mean = mean, level = level, lambda = lambda,
...)
mean.se = y.pca$mean.se
eigen_values = y.pca$eigen_value
}
if (weight == TRUE) {
q <- beta * (1 - beta)^(0:(n - 1))
my <- apply(y$y, 1, weighted.mean, w = rev(q))
newy <- sweep(t(y$y), 2, my)
x = y$x
y2 = y$y
n = dim(y$y)[2]
wq = diag(rev(q))
newy2 = wq %*% newy
dummy = svd(newy2)
load = dummy$v[, 1:(order)]
sco = newy %*% load
yy <- matrix(NA, nrow = ngrid, ncol = n)
xx = seq(min(x), max(x), l = ngrid)
for (i in 1:n) {
miss <- is.na(newy[i, ])
yy[, i] <- spline(x[!miss], t(newy)[!miss, i], n = ngrid)$y
}
V <- rep(1, ngrid) * t(yy)
dummy <- La.svd(V)
eigen_values = dummy$d^2
varprop = (dummy$d[1:order])^2/(sum((dummy$d)^2))
mean.se = approx(xx, sqrt(apply(yy, 1, var)/n), xout = x)$y
}
ytsp <- tsp(y$time)
if (weight == TRUE) {
colmeanrm = matrix(colMeans(sco), dim(sco)[2], 1)
scomeanrm = sweep(sco, 2, colmeanrm)
coeff <- ts(cbind(rep(1, dim(y$y)[2]), scomeanrm), start = ytsp[1],
frequency = ytsp[3])
}
else {
y.coeff <- matrix(NA, nrow=ncol(y$y), ncol=ncol(y.pca$coeff))
y.coeff[yynotNA, ] <- y.pca$coeff
colnames(y.coeff) <- colnames(y.pca$coeff)
coeff <- ts(y.coeff, start = ytsp[1], frequency = ytsp[3])
}
if (weight == TRUE) {
my = my + load %*% colmeanrm
basis <- cbind(my, load)
if(order == 0)
{
colnames(basis) = "mean"
colnames(coeff) = "mean"
}
else
{
colnames(basis) = c("mean", paste("phi", 1:order, sep = ""))
colnames(coeff) = c("mean", paste("beta", 1:order, sep = ""))
}
fits <- fts(y$x, sweep(load %*% t(scomeanrm), 1, my,
"+"), start = ytsp[1], frequency = ytsp[3], xname = y$xname,
yname = paste("Fitted", y$yname))
}
else {
basis <- y.pca$basis
basis_obj = basis %*% t(coeff)
colnames(basis_obj) = colnames(y$y)
fits <- fts(y$x, basis_obj, start = ytsp[1],
frequency = ytsp[3], xname = y$xname, yname = paste("Fitted",
y$yname))
fits$x = y$x
if(order == 0)
{
colnames(basis) = "mean"
colnames(coeff) = "mean"
}
else
{
if(mean)
{
colnames(basis) = c("mean", paste("phi", 1:order, sep = ""))
colnames(coeff) = c("mean", paste("beta", 1:order, sep = ""))
}
else
{
colnames(basis) = c(paste("phi", 1:order, sep = ""))
colnames(coeff) = c(paste("beta", 1:order, sep = ""))
}
}
}
# rownames(basis) <- paste(y$x)
res <- fts(y$x, y$y - fits$y, start = ytsp[1], frequency = ytsp[3],
xname = y$xname, yname = paste("Residuals", y$yname))
if (weight == FALSE) {
out <- list(x1 = as.numeric(colnames(y$y)), y1 = as.numeric(rownames(y$y)),
y = fts(y$x, y$y, xname = y$xname, yname = y$yname),
basis = basis, coeff = coeff, fitted = fits, residuals = res,
varprop = y.pca$varprop, eigen_values = eigen_values,
wt = ts(y.pca$weights, start = ytsp[1],
frequency = ytsp[3]), v = ts(y.pca$v, start = ytsp[1],
frequency = ytsp[3]), basis2 = y.pca$basis2,
coeff2 = y.pca$coeff2, mean.se = mean.se, call = match.call())
}
else {
out <- list(x1 = as.numeric(colnames(y$y)), y1 = as.numeric(rownames(y$y)),
y = fts(y$x, y2, xname = y$xname, yname = y$yname),
basis = basis, coeff = coeff, fitted = fits, residuals = res,
varprop = varprop, eigen_values = eigen_values, wt = rev(q), mean.se = mean.se,
call = match.call())
}
return(structure(out, class = c("ftsm", "fm")))
}