Nothing
#' Get initial parameter estimates for estimation
#'
#' @param y an object created from stsm_detect_frequency
#' @param freq Frequency of the data
#' @param trend Trend specification ("random-walk", "random-walk-drift", "double-random-walk", "random-walk2").
#' @param cycle The period for the longer-term cycle
#' @param decomp Decomposition model ("tend-cycle-seasonal", "trend-seasonal", "trend-cycle", "trend-noise")
#' @param seasons The seasonal lengths to split the seasonality into
#' @param prior A data table created by stsm_prior
#' @param sig_level Significance level for statistical tests
#' @param arma Named vector with values for p and q corresponding to the ARMA(p,q) specification if
#' @param exo Matrix of exogenous variables. Can be used to specify regression effects or other seasonal effects like holidays, etc.
#' @param state_eqns Character vector of equations to apply exo_state to the unobserved components. If left as the default, then all variables in
#' exo_state will be applied to all the unobserved components. The equations should look like:
#' "trend ~ var - 1", "drift ~ var - 1", "cycle ~ var - 1", "seasonal ~ var - 1".
#' If only some equations are specified, it will be assumed that the exogenous data will be applied to only those specified equations.
#' @param interpolate Character string giving frequency to interpolate to: i.e. "quarterly", "monthly", "weekly", "daily"
#' cycle is set to 'arma'. If NA, then will auto-select the order.
#' @param interpolate_method Character string giving the interpolation method:
#' @import data.table
#' @return named vector containing the initial parameter estimates for estimation
stsm_init_pars = function(y, freq, trend, cycle, decomp = "", seasons = NULL, prior = NULL, sig_level = 0.01,
arma = c(p = NA, q = NA), exo = NULL, state_eqns = NULL,
interpolate = NA, interpolate_method = NA){
if(is.null(prior)){
prior = stsm_prior(y, freq, decomp, seasons)
}else{
prior = copy(prior)
}
#Trend and drift parameters
div = max(c(1, sum(sapply(c("cycle", "seasonal", "trend"), function(x){grepl(x, decomp)}))))
div = ifelse(trend == "random-walk-drift", div + 1, div)
r = prior$remainder/div
if(trend == "random-walk"){
par = c(sig_t = stats::sd(diff(prior$trend + r), na.rm = TRUE))
}else if(trend == "random-walk2"){
par = c(sig_t = stats::sd(diff(diff(prior$trend + r)), na.rm = TRUE))
}else if(trend %in% c("random-walk-drift", "double-random-walk")){
if(trend == "random-walk-drift"){
arima = tryCatch(forecast::Arima(prior$drift + r, order = c(1, 0, 0), include.constant = TRUE),
error = function(err){NULL})
if(!is.null(arima)){
arima = c(arima$coef, sigma = sqrt(arima$sigma2))
#Keep initial parameters away from the feasible borders
arima["ar1"] = ifelse(arima["ar1"] > 0.99, 0.75,
ifelse(arima["ar1"] < 0, 0.25, arima["ar1"]))
}else{
arima = c(sigma = stats::sd(prior$drift + r, na.rm = TRUE),
intercept = mean(prior$drift, na.rm = TRUE)*(1 - 0.75), ar1 = 0.75)
}
par = c(sig_t = unname(sqrt(abs(stats::var(diff(prior$trend + r), na.rm = TRUE) - arima["sigma"]^2/(1 - arima["ar1"]^2)))),
sig_d = unname(arima["sigma"]), d = unname(arima["intercept"]), phi_d = unname(arima["ar1"]))
rm(arima)
}else if(trend == "double-random-walk"){
par = c(sig_t = stats::sd(diff(prior$trend + r), na.rm = TRUE)/sqrt(2),
sig_d = stats::sd(diff(prior$trend + r), na.rm = TRUE)/sqrt(2))
}
}
#Cycle parameters
if(grepl("cycle", decomp) & ifelse(!is.null(cycle) & length(cycle) > 0, cycle != "arma", FALSE)){
arima = tryCatch(forecast::Arima(stats::ts(prior$cycle + r, frequency = freq), order = c(2, 0, 1)),
error = function(err){NULL})
if(is.null(cycle)){
cycle = freq*5
}
#Set phi_c to be related to the speed of convergence of cycle towards steady state: i.e. maximum eigenvalue of ARMA(2,1)
if(!is.null(arima)){
Fm = rbind(arima$coef[c("ar1", "ar2", "ma1")],
c(1, 0, 0),
c(0, 0, 0))
phi_c = max(Mod(eigen(Fm)$values))
phi_c = max(c(0.01, min(c(phi_c, 1))))
sig_c = sqrt(arima$sigma2)
}else{
phi_c = 0.95
sig_c = stats::sd(diff(prior$cycle + r), na.rm = TRUE)
}
par = c(par, phi_c = phi_c, lambda = 2*pi/cycle, sig_c = sig_c)
}else if(grepl("cycle", decomp) & ifelse(!is.null(cycle) & length(cycle) > 0, cycle == "arma", FALSE)){
if(any(sapply(arma, is.na))){
if(!is.na(interpolate)){
arima = tryCatch(forecast::auto.arima(stats::ts(prior[!is.na(y), ]$cycle + r[!is.na(y)], frequency = freq),
max.d = 0, seasonal = FALSE, stationary = TRUE),
error = function(err){NULL})
if(!is.null(arima)){
arima = tryCatch(forecast::Arima(stats::ts(prior$cycle + r, frequency = freq), include.mean = any(grepl("mean|intercept", names(arima$coef))),
order = c(sum(grepl("ar\\d+", names(arima$coef))), 0, sum(grepl("ma\\d+", names(arima$coef)))), seasonal = c(0, 0, 0)),
error = function(err){NULL})
}
}else{
arima = tryCatch(forecast::auto.arima(stats::ts(prior$cycle + r, frequency = freq),
max.d = 0, seasonal = FALSE, stationary = TRUE),
error = function(err){NULL})
}
}else{
arima = tryCatch(forecast::Arima(stats::ts(prior$cycle + r, frequency = freq),
order = c(arma["p"], 0, arma["q"]), seasonal = c(0, 0, 0)),
error = function(err){NULL})
}
if(!is.null(arima) & length(arima$coef) > 0){
par = c(par, arima$coef[grepl("ar\\d+|ma\\d+", names(arima$coef))], sig_c = sqrt(arima$sigma2))
names(par)[grepl("ar\\d+", names(par))] = gsub("ar", "phi_c.", names(par)[grepl("ar\\d+", names(par))])
names(par)[grepl("ma\\d+", names(par))] = gsub("ma", "theta_c.", names(par)[grepl("ma\\d+", names(par))])
}else{
par = c(par, phi_c.1 = 1.25, phi_c.2 = -0.3, theta_c.1 = 0.75, sig_c = stats::sd(diff(prior$cycle + r)))
}
}
#Seasonal parameters
if(grepl("seasonal", decomp)){
par = c(par, sig_s = unname(rep(stats::sd(diff(prior$seasonal + r), na.rm = TRUE)/sqrt(length(seasons)), length(seasons))))
names(par)[grepl("sig_s", names(par))] = paste0("sig_s", seasons)
}
par = c(par, sig_e = stats::sd(r, na.rm = TRUE))
if(any(par[grepl("sig_", names(par))] <= 0)){
par[grepl("sig_", names(par))] = ifelse(stats::sd(prior$remainder, na.rm = TRUE) > 0, stats::sd(prior$remainder, na.rm = TRUE), 0.01)
}
#Exogenous parameters
if(is.null(exo)){
par = c(par, betaO_Xo = NA, betaS_Xs = NA)
}else{
if(sum(grepl("obs\\.", colnames(exo))) > 0){
lm_obs = stats::lm(y ~ . - 1, data = data.frame(cbind(y,
prior[, colnames(prior) %in% c("trend", "cycle", "seasonal"), with = FALSE],
exo[, grepl("obs\\.", colnames(exo)), with = FALSE])))
coef_obs = stats::coef(lm_obs)
coef_obs = coef_obs[grepl("obs\\.", names(coef_obs))]
names(coef_obs) = gsub("obs\\.", "", names(coef_obs))
names(coef_obs) = paste0("betaO_", names(coef_obs))
par = c(par, coef_obs)
}else{
par = c(par, betaO_Xo = NA)
}
if(sum(grepl("state\\.", colnames(exo))) > 0){
ssm = stsm_ssm(par, y, decomp, trend, init = NULL, prior = prior, freq = freq, seasons = seasons,
interpolate = interpolate, interpolate_method = interpolate_method)
Fm = ssm[["Fm"]]
colnames(Fm)[grepl("^Tt_\\d+$", colnames(Fm))] =
gsub("^Tt_", "trend.l", colnames(Fm)[grepl("^Tt_\\d+$", colnames(Fm))])
colnames(Fm)[grepl("^Dt_\\d+$", colnames(Fm))] =
gsub("^Dt_", "drift.l", colnames(Fm)[grepl("^Dt_\\d+$", colnames(Fm))])
colnames(Fm)[grepl("^Ct_\\d+$", colnames(Fm))] =
gsub("^Ct_", "cycle.l", colnames(Fm)[grepl("^Ct_\\d+$", colnames(Fm))])
colnames(Fm)[grepl("^St\\d+", colnames(Fm))] =
gsub("\\_", "\\.l", gsub("St", "seasonal", colnames(Fm)[grepl("^St\\d+", colnames(Fm))]))
lag = max(unlist(lapply(strsplit(colnames(Fm), "\\.l"), function(x){as.numeric(x[2])})), na.rm = TRUE)
cols = copy(colnames(prior))
for(i in 1:lag){
prior[, paste0(cols, ".l", i) := lapply(.SD, function(x){
shift(x, type = "lag", n = i)
}), .SDcols = c(cols)]
}
if(is.null(state_eqns)){
state_eqns = c("trend ~ . - 1", "drift ~ . - 1", "cycle ~ . - 1", "seasonal ~ . - 1")
}else{
for(i in 1:length(state_eqns)){
names(state_eqns)[i] = trimws(strsplit(state_eqns[i], "~")[[1]][1])
}
}
rownames = rownames(Fm)[grepl("_0", rownames(Fm)) & !grepl("Sts|Cts", rownames(Fm))]
mat = lapply(rownames, function(x){
var = ifelse(x == "Tt_0", "trend",
ifelse(x == "Dt_0", "drift",
ifelse(x == "Ct_0", "cycle",
ifelse(grepl("St", x), gsub("_0", "", gsub("St", "seasonal", x)), NA))))
var2 = gsub("[[:digit:]]|\\.", "", var)
eqn = gsub("\\.", paste(colnames(exo)[grepl("state\\.", colnames(exo))], collapse = " + "), state_eqns[var2])
if(is.na(eqn)){
return(data.table(NULL))
}else{
if(!grepl("\\-1|\\- 1", eqn)){
eqn = paste(eqn, "- 1")
}
eqn_split = trimws(strsplit(eqn, "~|\\+|\\-")[[1]])
eqn = paste0(eqn_split[1], " ~ ", paste(unlist(lapply(eqn_split[2:length(eqn_split)], function(x){
if(trimws(x) != "1"){
x = paste0("state.", trimws(x))
}
return(x)
})), collapse = " + "))
eqn = gsub("\\+1|\\+ 1", "- 1", eqn)
}
eqn2 = c(names(which(Fm[x, ] != 0 & !grepl("^Cts|^Sts", colnames(Fm)))))
eqn = paste(trimws(strsplit(eqn, "~")[[1]][1]), "~",
paste(eqn2, collapse = " + "), "+",
trimws(strsplit(eqn, "~")[[1]][2]))
dt = as.data.frame(cbind(prior, exo))
lm = stats::lm(eqn, data = dt)
return(list(name = x, table = as.data.table(t(stats::coef(lm)[grepl("state\\.", names(stats::coef(lm)))]))))
})
mat = rbindlist(lapply(mat, function(x){x[["table"]]}), use.names = TRUE, fill = TRUE)[, "rn" := unlist(lapply(mat, function(x){x[["name"]]}))]
mat = as.matrix(merge(data.table(rn = rownames(Fm)), mat, by = "rn", all = TRUE, sort = FALSE), rownames = "rn")
colnames(mat) = gsub("state\\.", "", colnames(mat))
if("Cts_0" %in% rownames(mat)){
mat["Cts_0", ] = mat["Ct_0", ]
}
mat[grepl("^Sts\\d+", rownames(mat)), ] = mat[grepl("^St\\d+", rownames(mat)), ]
cols = unique(gsub("state\\.", "", colnames(exo)[grepl("state\\.", colnames(exo))]))
par = c(par, unname(c(mat[, colnames(mat) %in% cols])))
names(par)[names(par) == ""] = unlist(lapply(colnames(mat)[colnames(mat) %in% cols], function(x){
paste0("betaS_", rownames(mat), ".", x)
}))
}else{
par = c(par, betaS_Xs = NA)
}
}
return(par)
}
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.