dynamic_FLR <- function (dat, newdata, holdoutdata, order_k_percent = 0.9, order_m_percent = 0.9,
pcd_method = c("classical", "M"), robust_lambda = 2.33, bootrep = 100,
pointfore, level = 80)
{
newdata = matrix(newdata, nrow = 1)
holdoutdata = matrix(holdoutdata, nrow = 1)
pcd_method = match.arg(pcd_method)
if (ncol(newdata) == 1) {
data_first = t(scale(dat[1, ], scale = FALSE))
data_last = t(scale(t(dat[(ncol(newdata) + 1):nrow(dat),
]), center = TRUE, scale = FALSE))
}
if (ncol(holdoutdata) == 1) {
data_first = t(scale(t(dat[1:ncol(newdata), ]), center = TRUE,
scale = FALSE))
data_last = t(scale(dat[nrow(dat), ], center = TRUE,
scale = FALSE))
}
if (ncol(newdata) > 1 & ncol(holdoutdata) > 1) {
data_first = t(scale(t(dat[1:ncol(newdata), ]), center = TRUE,
scale = FALSE))
data_last = t(scale(t(dat[(ncol(newdata) + 1):nrow(dat),
]), center = TRUE, scale = FALSE))
}
cross_cov = data_last %*% t(data_first)
order_k = head(which(round(cumsum(eigen(crossprod(t(data_first)))$value)/sum((eigen(crossprod(t(data_first)))$value)),
2) >= order_k_percent), 1)
if (pcd_method == "classical") {
phi_k = matrix(eigen(crossprod(t(data_first)))$vectors[, 1:order_k], ncol = order_k)
xi_k = t(data_first) %*% phi_k
reconstruct_k = phi_k %*% t(xi_k)
resi_k = data_first - reconstruct_k
lambda_k = eigen(crossprod(t(data_first)))$values[1:order_k]
if (pointfore == FALSE) {
n_xi = nrow(xi_k)
p_xi = ncol(xi_k)
if (requireNamespace("meboot", quietly = TRUE)) {
out_xi <- meboot::meboot(ts(as.numeric(xi_k),
start = 1, end = n_xi * p_xi), reps = bootrep)
}
else {
stop("Please install meboot")
}
n_resi = nrow(t(resi_k))
p_resi = ncol(t(resi_k))
out_resi <- meboot::meboot(ts(as.numeric(t(resi_k)),
start = 1, end = n_resi * p_resi), reps = bootrep)
data_first_boot = array(NA, dim = c(nrow(data_first),
ncol(data_first), bootrep))
phi_k_boot = array(NA, dim = c(nrow(phi_k), ncol(phi_k),
bootrep))
lambda_k_boot = matrix(NA, order_k, bootrep)
for (ik in 1:bootrep) {
xi_k_boot = matrix(out_xi$ensemble[, ik], n_xi,
p_xi)
resi_k_boot = t(matrix(out_resi$ensemble[, ik],
n_resi, p_resi))
if (!any(is.finite(resi_k_boot))) {
data_first_boot[, , ik] = phi_k %*% t(xi_k_boot)
}
else {
data_first_boot[, , ik] = phi_k %*% t(xi_k_boot) +
resi_k_boot
}
phi_k_boot[, , ik] = matrix(eigen(crossprod(t(data_first_boot[, , ik])))$vectors[, 1:order_k], ncol = order_k)
lambda_k_boot[, ik] = eigen(crossprod(t(data_first_boot[,, ik])))$values[1:order_k]
}
}
}
if (pcd_method == "M") {
if (length(newdata) == 1) {
phi_k = matrix(eigen(crossprod(t(data_first)))$vectors[,
1:order_k], ncol = order_k)
}
else {
phi_k = matrix(ftsm(fts(1:length(newdata), data_first),
order = order_k, method = "M", lambda = robust_lambda,
mean = FALSE)$basis[, 1:order_k], ncol = order_k)
}
xi_k = t(data_first) %*% phi_k
reconstruct_k = phi_k %*% t(xi_k)
resi_k = data_first - reconstruct_k
lambda_k = eigen(crossprod(t(data_first)))$values[1:order_k]
if (pointfore == FALSE) {
n_xi = nrow(xi_k)
p_xi = ncol(xi_k)
out_xi <- meboot::meboot(ts(as.numeric(xi_k), start = 1,
end = n_xi * p_xi), reps = bootrep)
n_resi = nrow(t(resi_k))
p_resi = ncol(t(resi_k))
out_resi <- meboot::meboot(ts(as.numeric(t(resi_k)),
start = 1, end = n_resi * p_resi), reps = bootrep)
data_first_boot = array(NA, dim = c(nrow(data_first),
ncol(data_first), bootrep))
phi_k_boot = array(NA, dim = c(nrow(phi_k), ncol(phi_k),
bootrep))
lambda_k_boot = matrix(NA, order_k, bootrep)
for (ik in 1:bootrep) {
xi_k_boot = matrix(out_xi$ensemble[, ik], n_xi,
p_xi)
resi_k_boot = t(matrix(out_resi$ensemble[, ik],
n_resi, p_resi))
if (!any(is.finite(resi_k_boot))) {
data_first_boot[, , ik] = phi_k %*% t(xi_k_boot)
}
else {
data_first_boot[, , ik] = phi_k %*% t(xi_k_boot) +
resi_k_boot
}
phi_k_boot[, , ik] = matrix(eigen(crossprod(t(data_first_boot[,
, ik])))$vectors[, 1:order_k], ncol = order_k)
lambda_k_boot[, ik] = eigen(crossprod(t(data_first_boot[,
, ik])))$values[1:order_k]
}
}
}
order_m = head(which(round(cumsum((eigen(crossprod(t(data_last)))$value)/sum((eigen(crossprod(t(data_last)))$value))),
2) >= order_m_percent), 1)
if (pcd_method == "classical") {
psi_m = as.matrix(eigen(crossprod(t(data_last)))$vectors[,
1:order_m])
eta_m = t(data_last) %*% psi_m
reconstruct_m = psi_m %*% t(eta_m)
resi_m = data_last - reconstruct_m
if (pointfore == FALSE) {
n_eta = nrow(eta_m)
p_eta = ncol(eta_m)
out_eta = meboot::meboot(ts(as.numeric(eta_m), start = 1,
end = n_eta * p_eta), reps = bootrep)
n_resi = nrow(t(resi_m))
p_resi = ncol(t(resi_m))
out_resi = meboot::meboot(ts(as.numeric(t(resi_m)),
start = 1, end = n_resi * p_resi), reps = bootrep)
data_last_boot = array(NA, dim = c(nrow(data_last),
ncol(data_last), bootrep))
psi_m_boot = array(NA, dim = c(nrow(psi_m), ncol(psi_m),
bootrep))
for (ik in 1:bootrep) {
eta_m_boot = matrix(out_eta$ensemble[, ik], n_eta,
p_eta)
resi_m_boot = t(matrix(out_resi$ensemble[, ik],
n_resi, p_resi))
if (!any(is.finite(resi_m_boot))) {
data_last_boot[, , ik] = psi_m %*% t(eta_m_boot)
}
else {
data_last_boot[, , ik] = psi_m %*% t(eta_m_boot) +
resi_m_boot
}
if (nrow(data_last) == 1) {
psi_m_boot[, , ik] = matrix(eigen(crossprod(data_last_boot[,
, ik]))$vectors[, 1:order_m], ncol = order_m)
}
else {
psi_m_boot[, , ik] = matrix(eigen(crossprod(t(data_last_boot[,
, ik])))$vectors[, 1:order_m], ncol = order_m)
}
}
}
}
if (pcd_method == "M") {
if (length(holdoutdata) == 1) {
psi_m = as.matrix(eigen(crossprod(t(data_last)))$vectors[,
1:order_m])
}
else {
psi_m = as.matrix(ftsm(fts(1:length(holdoutdata),
data_last), order = order_m, method = "M", lambda = robust_lambda,
mean = FALSE)$basis[, 1:order_m])
}
eta_m = t(data_last) %*% psi_m
reconstruct_m = psi_m %*% t(eta_m)
resi_m = data_last - reconstruct_m
if (pointfore == FALSE) {
n_eta = nrow(eta_m)
p_eta = ncol(eta_m)
out_eta = meboot::meboot(ts(as.numeric(eta_m), start = 1,
end = n_eta * p_eta), reps = bootrep)
n_resi = nrow(t(resi_m))
p_resi = ncol(t(resi_m))
out_resi = meboot::meboot(ts(as.numeric(t(resi_m)),
start = 1, end = n_resi * p_resi), reps = bootrep)
data_last_boot = array(NA, dim = c(nrow(data_last),
ncol(data_last), bootrep))
psi_m_boot = array(NA, dim = c(nrow(psi_m), ncol(psi_m),
bootrep))
for (ik in 1:bootrep) {
eta_m_boot = matrix(out_eta$ensemble[, ik], n_eta,
p_eta)
resi_m_boot = t(matrix(out_resi$ensemble[, ik],
n_resi, p_resi))
if (!any(is.finite(resi_m_boot))) {
data_last_boot[, , ik] = psi_m %*% t(eta_m_boot)
}
else {
data_last_boot[, , ik] = psi_m %*% t(eta_m_boot) +
resi_m_boot
}
if (nrow(data_last) == 1) {
psi_m_boot[, , ik] = matrix(eigen(crossprod(data_last_boot[,
, ik]))$vectors[, 1:order_m], ncol = order_m)
}
else {
psi_m_boot[, , ik] = matrix(eigen(crossprod(t(data_last_boot[,
, ik])))$vectors[, 1:order_m], ncol = order_m)
}
}
}
}
if (pointfore == TRUE) {
lambda_km = matrix(NA, order_k, order_m)
for (j in 1:order_k) {
for (i in 1:order_m) {
lambda_km[j, i] = matrix(psi_m[, i], nrow = 1) %*%
cross_cov %*% matrix(phi_k[, j], nrow = length(newdata))
}
}
beta = array(NA, dim = c(length(newdata), length(holdoutdata),
order_m, order_k))
for (j in 1:order_k) {
for (i in 1:order_m) {
beta[, , i, j] = lambda_km[j, i]/lambda_k[j] *
matrix(phi_k[, j], ncol = 1) %*% matrix(psi_m[,
i], nrow = 1)
}
}
finalbeta = apply(beta, c(1, 2), sum)
if (ncol(newdata) == 1) {
update_forecast = rowMeans(dat[(ncol(newdata) + 1):nrow(dat),
]) + matrix(newdata - mean(dat[1:ncol(newdata),
]), nrow = 1) %*% finalbeta
}
if (ncol(holdoutdata) == 1) {
update_forecast = mean(dat[(ncol(newdata) + 1):nrow(dat),
]) + matrix(newdata - rowMeans(dat[1:ncol(newdata),
]), nrow = 1) %*% finalbeta
}
if (ncol(newdata) > 1 & ncol(holdoutdata) > 1) {
update_forecast = rowMeans(dat[(ncol(newdata) + 1):nrow(dat),
]) + matrix(newdata - rowMeans(dat[1:ncol(newdata),
]), nrow = 1) %*% finalbeta
}
err = matrix(c(error(forecast = update_forecast, true = holdoutdata,
method = "mae"), error(forecast = update_forecast,
true = holdoutdata, method = "mse")), nrow = 1)
colnames(err) = c("MAFE", "MSFE")
return(list(update_forecast = update_forecast, holdoutdata = holdoutdata,
err = err, order_k = order_k, order_m = order_m))
}
else {
lambda_km_boot = array(NA, dim = c(order_k, order_m, bootrep))
for (ik in 1:bootrep) {
for (j in 1:order_k) {
for (i in 1:order_m) {
cross_cov = data_last_boot[, , ik] %*% t(data_first_boot[,
, ik])
lambda_km_boot[j, i, ik] = matrix(psi_m_boot[,
i, ik], nrow = 1) %*% cross_cov %*% matrix(phi_k_boot[,
j, ik], nrow = length(newdata))
}
}
}
beta_boot = array(NA, dim = c(length(newdata), length(holdoutdata),
order_m, order_k, bootrep))
for (ik in 1:bootrep) {
for (j in 1:order_k) {
for (i in 1:order_m) {
beta_boot[, , i, j, ik] = lambda_km_boot[j,
i, ik]/lambda_k_boot[j, ik] * matrix(phi_k_boot[,
j, ik], ncol = 1) %*% matrix(psi_m_boot[,
i, ik], nrow = 1)
}
}
}
update_forecast = matrix(NA, nrow(dat) - length(newdata), bootrep)
for (ik in 1:bootrep) {
if (ncol(holdoutdata) == 1)
{
finalbeta_boot = beta_boot[, , , , ik]
update_forecast[, ik] = mean(dat[(ncol(newdata) + 1):nrow(dat), ]) + matrix(newdata - rowMeans(dat[1:ncol(newdata),]), nrow = 1) %*% finalbeta_boot
}
if (ncol(newdata) == 1) {
finalbeta_boot = apply(beta_boot[, , , , ik], c(1, 2), sum)
update_forecast[, ik] = rowMeans(dat[(ncol(newdata) +
1):nrow(dat), ]) + matrix(newdata - mean(dat[1:ncol(newdata),
]), nrow = 1) %*% finalbeta_boot
}
if (ncol(newdata) > 1 & ncol(holdoutdata) > 1) {
finalbeta_boot = apply(beta_boot[, , , , ik], c(1, 2), sum)
update_forecast[, ik] = rowMeans(dat[(ncol(newdata) +
1):nrow(dat), ]) + matrix(newdata - rowMeans(dat[1:ncol(newdata),
]), nrow = 1) %*% finalbeta_boot
}
}
ik = nrow(dat) - length(holdoutdata)
err = matrix(NA, length(holdoutdata), ncol(dat) - 2)
for (j in 2:(ncol(dat) - 1)) {
holdout = dat[(ik + 1):nrow(dat), (j + 1)]
dum = dynamic_FLR(dat = dat[, 1:j], newdata = dat[1:ik,
(j + 1)], holdoutdata = holdout, order_k_percent = order_k_percent,
order_m_percent = order_m_percent, pointfore = TRUE)
err[, j - 1] = holdout - dum$update_forecast
}
err_boot = err[, sample(1:(ncol(dat) - 2), bootrep, replace = TRUE)]
update_comb = update_forecast + err_boot
update_comb_lb_ub = apply(update_comb, 1, quantile, c((100 -
level)/200, (100 + level)/200))
return(list(update_comb = update_comb, update_comb_lb_ub = update_comb_lb_ub,
err_boot = err_boot, update_forecast = update_forecast))
}
}