## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup, warning=FALSE,message=FALSE---------------------------------------
library(tsissm)
library(tsaux)
library(xts)
library(data.table)
library(zoo)

## -----------------------------------------------------------------------------
data("us_retail_sales")
y <- as.xts(us_retail_sales)
train <- y["/2023"]
test <- y["2024/"]
lambda_pre_estimate <- box_cox(lambda = NA)$transform(train) |> attr("lambda")
xreg <- auto_regressors(train, frequency = 12, lambda = lambda_pre_estimate, sampling = "months", h = 14, method = "sequential", 
                        check.rank = TRUE, discard.cval = 3.5, maxit.iloop = 10, maxit.oloop = 10, 
                        forc_dates = index(test))
spec <- issm_modelspec(train, slope = TRUE, seasonal = TRUE, seasonal_frequency = 12,
                       seasonal_harmonics = 5, ar = 1, ma = 0, lambda = lambda_pre_estimate, 
                       xreg = xreg$xreg[index(train),])
spec$parmatrix[grepl("^kappa", parameters), initial := xreg$init]
mod <- spec |> estimate()

## -----------------------------------------------------------------------------
test_xreg <- xreg$xreg[index(test),]
predict_list <- vector(mode = "list", length = 14)
predict_list[[1]] <- predict(mod, h = 1, newxreg = test_xreg[1,])
new_mod <- mod
for (i in 2:14) {
    new_mod <-  tsfilter(new_mod, y = test[i - 1], newxreg = test_xreg[i - 1, ])
    predict_list[[i]] <- predict(new_mod, h = 1, newxreg = test_xreg[i,])
}
forecast_distribution <- do.call(cbind, lapply(predict_list, function(x) x$distribution))
forecast_mean <- do.call(rbind, lapply(predict_list, function(x) x$analytic_mean))
class(forecast_distribution) <- "tsmodel.distribution"
test_transformed <- cbind(fitted(new_mod)[index(test)], forecast_mean)
colnames(test_transformed) <- c("Filtered","Forecast")
head(test_transformed)

## -----------------------------------------------------------------------------
p <- tsmoments(mod, h = 1, newxreg = test_xreg[1,], transform = FALSE)
tmp_mod <-  tsfilter(mod, y = test[1], newxreg = test_xreg[1, ])
fitted_value <- tail(tmp_mod$spec$target$y, 1) - tail(tmp_mod$model$error,1)
all.equal(p$mean, fitted_value)

## ----fig.width=6,fig.height=4-------------------------------------------------
filtered_model <- fitted(new_mod)
original_model <- fitted(mod)
forecast_model <- forecast_mean
Z <- cbind(original_model, filtered_model, forecast_model)
matplot(index(tail(Z,50)), coredata(tail(Z,50)), type = c("l","l","p"), col = c("black","orange","steelblue"), pch = c(1,1,1), lty = c(1,2,1), ylab = "", xlab = "", lwd = c(1,1.5, 1.8))
legend("topleft", c("Fitted","Filtered","Forecast"), col = c("black","orange","steelblue"), lty = c(1,1,NA), lwd = c(1,1.5, 1.8), pch = c(NA,NA,1), bty = "n")