## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 3.4,
                      message = FALSE, warning = FALSE)
options(digits = 4)

## ----setup--------------------------------------------------------------------
library(RSDC)

## ----data---------------------------------------------------------------------
set.seed(1)
n  <- 300L
X  <- cbind(intercept = 1, z = as.numeric(scale(seq_len(n))))   # include an intercept
beta <- rbind(c(1.5, -0.6),     # regime 1 (logistic stay-prob coefficients)
              c(1.0,  0.2))     # regime 2
rho  <- rbind(0.2, 0.8)         # K = 2 -> one correlation per regime
Sig  <- array(c(matrix(c(1, 0.2, 0.2, 1), 2),
                matrix(c(1, 0.8, 0.8, 1), 2)), dim = c(2, 2, 2))
sim  <- rsdc_simulate(n = n, X = X, beta = beta, mu = matrix(0, 2, 2),
                      sigma = Sig, N = 2, seed = 42)
y <- scale(sim$observations)    # treated as unit-variance standardized residuals

ctrl <- list(itermax = 30, NP = 40, seed = 1)   # small budget for the vignette
fit <- rsdc_estimate("tvtp", residuals = y, N = 2, X = X, control = ctrl)
class(fit)
fit

## ----methods------------------------------------------------------------------
logLik(fit)
c(AIC = AIC(fit), BIC = BIC(fit), nobs = nobs(fit))
coef(fit)

## ----summary------------------------------------------------------------------
summary(fit)
# Standard errors are unavailable when the optimum sits on the feasible
# boundary, so guard the interval call rather than assume it succeeds.
if (!is.null(fit$vcov)) head(confint(fit))

## ----compare------------------------------------------------------------------
fit_const <- rsdc_estimate("const", residuals = y, control = ctrl)
fit_noX   <- rsdc_estimate("noX",   residuals = y, N = 2, control = ctrl)
data.frame(
  model = c("const", "noX(N=2)", "tvtp(N=2)"),
  AIC   = c(AIC(fit_const), AIC(fit_noX), AIC(fit)),
  BIC   = c(BIC(fit_const), BIC(fit_noX), BIC(fit))
)

## ----predict------------------------------------------------------------------
sigma_mat <- matrix(1, n, 2, dimnames = list(NULL, c("a", "b")))   # unit vols here
fc <- predict(fit, residuals = y, sigma_matrix = sigma_mat,
              value_cols = c("a", "b"), X = X)
str(fc$predicted_correlations)

## ----simulate-----------------------------------------------------------------
sim2 <- simulate(fit, X = X, seed = 7)
table(sim2$states)

## ----oos----------------------------------------------------------------------
fc_oos <- predict(fit, residuals = y, sigma_matrix = sigma_mat,
                  value_cols = c("a", "b"), X = X, out_of_sample = TRUE)
nrow(fc_oos$predicted_correlations)   # length of the 30% hold-out
fc_oos$BIC                            # out-of-sample predictive score

## ----multistart---------------------------------------------------------------
fit_ms <- rsdc_estimate("tvtp", residuals = y, N = 2, X = X,
                        control = c(ctrl, list(n_starts = 3)))
fit_ms$start_logliks            # one log-likelihood per start

## ----bands--------------------------------------------------------------------
# Bands are drawn from the estimated sampling distribution, so they need a
# covariance; skip gracefully if this quick fit did not produce one.
if (!is.null(fit$vcov)) {
  bands <- rsdc_corr_bands(fit, B = 100, seed = 1)
  head(bands[[1]])
}

## ----viterbi------------------------------------------------------------------
table(rsdc_viterbi(fit))

## ----ahead--------------------------------------------------------------------
fa <- rsdc_forecast_ahead(fit, horizon = 5)
fa$predicted_correlations

## ----broom--------------------------------------------------------------------
generics::tidy(fit)
generics::glance(fit)

## ----autoplot, eval = requireNamespace("ggplot2", quietly = TRUE), fig.alt = "Stacked smoothed regime probabilities over time"----
ggplot2::autoplot(fit)          # stacked smoothed regime probabilities

## ----bign---------------------------------------------------------------------
fit4 <- rsdc_estimate("noX", residuals = y, N = 4,
                      control = list(itermax = 20, NP = 60, seed = 1))
dim(fit4$transition_matrix)

## ----recovery-----------------------------------------------------------------
rbind(true = c(0.2, 0.8),
      estimated = sort(as.numeric(fit$correlations)))