https://github.com/cran/ftsa
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
dynamic_FLR.R
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)
#############
# FPCA (1st)
#############
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(, dim = c(nrow(data_first), ncol(data_first), bootrep))
phi_k_boot = array(, dim = c(nrow(phi_k), ncol(phi_k), bootrep))
lambda_k_boot = matrix(, 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(, dim = c(nrow(data_first), ncol(data_first), bootrep))
phi_k_boot = array(, dim = c(nrow(phi_k), ncol(phi_k), bootrep))
lambda_k_boot = matrix(, 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]
}
}
}
#############
# FPCA (2nd)
#############
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(, dim = c(nrow(data_last), ncol(data_last), bootrep))
psi_m_boot = array(, 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(, dim = c(nrow(data_last), ncol(data_last), bootrep))
psi_m_boot = array(, 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(, 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(, 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(, 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(, 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(, nrow(dat) - length(newdata), bootrep)
for(ik in 1:bootrep)
{
finalbeta_boot = apply(beta_boot[,,,,ik], c(1,2), sum)
if(ncol(newdata) == 1)
{
update_forecast[,ik] = rowMeans(dat[(ncol(newdata)+1):nrow(dat),]) + matrix(newdata - mean(dat[1:ncol(newdata),]), nrow=1)%*%finalbeta_boot
}
if(ncol(holdoutdata) == 1)
{
finalbeta_boot = as.matrix(rowSums(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 & ncol(holdoutdata) > 1)
{
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(, 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))
}
}