## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----eval=FALSE---------------------------------------------------------------
# install.packages("rbreak")

## -----------------------------------------------------------------------------
library(rbreak)

## -----------------------------------------------------------------------------
data(example)
y <- example$y
bigt <- length(y)
cat("Sample size:", bigt, "\n")

## ----fig.width=7, fig.height=4------------------------------------------------
plot(y, type = "l", 
     main = "Example Time Series Data",
     xlab = "Time", ylab = "Value", 
     col = "steelblue", lwd = 0.8)
grid()

## -----------------------------------------------------------------------------
# Load the example data
data(example)
y <- example$y
bigt <- length(y)

# Regressors (intercept only in this example)
z <- matrix(1, nrow = bigt, ncol = 1)

# Number of regressors
q <- 1

# Number of breaks under the alternative
m <- 3

## ----eval=FALSE---------------------------------------------------------------
# # Form A: delta = S*theta + ss
# # We have 4 regimes with 1 regressor each, so delta is 4x1
# # We want: delta_1 = theta_1, delta_2 = theta_2, delta_3 = theta_1, delta_4 = theta_2
# # This means we only have 2 free parameters (theta_1, theta_2)
# 
# S <- matrix(c(1, 0,    # delta_1 = theta_1
#               0, 1,    # delta_2 = theta_2
#               1, 0,    # delta_3 = theta_1
#               0, 1),   # delta_4 = theta_2
#             nrow = 4, ncol = 2, byrow = TRUE)
# ss <- matrix(0, nrow = 4, ncol = 1)

## ----eval=FALSE---------------------------------------------------------------
# result_formA <- mainp(
#   m = m, q = q, z = z, y = y,
#   trm = 0.10,          # 10% trimming
#   robust = 1,          # Use HAC standard errors
#   prewhit = 0,         # No prewhitening
#   hetvar = 1,          # Allow heteroskedastic variance
#   S = S, ss = ss,      # Form A restriction
#   R = NULL, rr = NULL, # Not used when forma=1
#   doestim = 1,         # Perform estimation
#   dotest = 1,          # Perform sup-F test
#   docv = 0,            # Critical value simulation
#   bigt = bigt,
#   forma = 1            # Use Form A restrictions
# )

## ----eval=FALSE---------------------------------------------------------------
# # Restriction matrix: R*delta = rr
# # delta = (delta_1, delta_2, delta_3, delta_4)
# # Restrictions: delta_1 = delta_3, delta_2 = delta_4
# R <- matrix(c(1, 0, -1, 0,   # delta_1 - delta_3 = 0
#               0, 1,  0, -1),  # delta_2 - delta_4 = 0
#             nrow = 2, byrow = TRUE)
# rr <- c(0, 0)

## ----eval=FALSE---------------------------------------------------------------
# result <- mainp(
#   m = m, q = q, z = z, y = y,
#   trm = 0.10,          # 10% trimming
#   robust = 1,          # Use HAC standard errors
#   prewhit = 0,         # No prewhitening
#   hetvar = 1,          # Allow heteroskedastic variance
#   R = R, rr = rr,      # Form B restriction
#   doestim = 1,         # Perform estimation
#   dotest = 1,          # Perform sup-F test
#   docv = 0,            # Critical value simulation
#   bigt = bigt,
#   formb = 1            # Use Form B restrictions
# )

## ----eval=FALSE---------------------------------------------------------------
# # Initial break dates (can be from unrestricted estimation or user-specified)
# # For the example data, we use manually specified initial dates
# T_init <- c(47, 64, 87)  # Initial break dates
# 
# # Run BRBP algorithm (Form B)
# result_brbp <- brbp(
#   y = y, z = z, q = q, m = m, bigt = bigt, trm = 0.10,
#   T_init = T_init,
#   B = 20,              # Number of bootstrap replications
#   R = R, rr = rr,      # Form B restrictions
#   verbose = FALSE      # Set to TRUE to see progress
# )
# 
# # Results
# breaks_brbp <- result_brbp$dx
# coefficients_brbp <- result_brbp$delta
# rssr_brbp <- result_brbp$rssr
# 
# # Track convergence
# all_breaks <- result_brbp$all_break_dates  # All break dates across iterations
# all_rssr <- result_brbp$all_rssr           # All RSSR values across iterations
# 
# # Plot convergence (if desired)
# plot(all_rssr, type = "l", xlab = "Iteration", ylab = "RSSR",
#      main = "BRBP Convergence (Form B)")

## ----eval=FALSE---------------------------------------------------------------
# # Form A restrictions (same as defined earlier)
# S <- matrix(c(1, 0, 1, 0, 0, 1, 0, 1), nrow = 4, ncol = 2, byrow = TRUE)
# ss <- matrix(0, nrow = 4, ncol = 1)
# 
# # Run BRBP algorithm (Form A)
# result_brbp_formA <- brbp(
#   y = y, z = z, q = q, m = m, bigt = bigt, trm = 0.10,
#   T_init = T_init,
#   B = 20,              # Number of bootstrap replications
#   S = S, ss = as.vector(ss),  # Form A restrictions
#   R = NULL, rr = NULL,        # Not used for Form A
#   verbose = FALSE
# )
# 
# # Results
# breaks_brbp_formA <- result_brbp_formA$dx
# coefficients_brbp_formA <- result_brbp_formA$delta
# rssr_brbp_formA <- result_brbp_formA$rssr
# 
# # Plot convergence
# plot(result_brbp_formA$all_rssr, type = "l", xlab = "Iteration", ylab = "RSSR",
#      main = "BRBP Convergence (Form A)")

## ----eval=FALSE---------------------------------------------------------------
# # Run RBRBP algorithm (Form B)
# result_rbrbp <- brbp_residual(
#   y = y, z = z, q = q, m = m, bigt = bigt, trm = 0.10,
#   T_init = T_init,
#   B = 20,              # Number of bootstrap replications
#   R = R, rr = rr,      # Form B restrictions
#   verbose = FALSE      # Set to TRUE to see progress
# )
# 
# # Results
# breaks_rbrbp <- result_rbrbp$dx
# coefficients_rbrbp <- result_rbrbp$delta
# rssr_rbrbp <- result_rbrbp$rssr
# 
# # Track convergence
# rall_breaks <- result_rbrbp$all_break_dates  # All break dates across iterations
# rall_rssr <- result_rbrbp$all_rssr           # All RSSR values across iterations
# 
# # Plot convergence (if desired)
# plot(rall_rssr, type = "l", xlab = "Iteration", ylab = "RSSR",
#      main = "RBRBP Convergence (Form B)")

## ----eval=FALSE---------------------------------------------------------------
# # Run RBRBP algorithm (Form A)
# result_rbrbp_formA <- brbp_residual(
#   y = y, z = z, q = q, m = m, bigt = bigt, trm = 0.10,
#   T_init = T_init,
#   B = 20,              # Number of bootstrap replications
#   S = S, ss = as.vector(ss),  # Form A restrictions
#   R = NULL, rr = NULL,        # Not used for Form A
#   verbose = FALSE
# )
# 
# # Results
# breaks_rbrbp_formA <- result_rbrbp_formA$dx
# coefficients_rbrbp_formA <- result_rbrbp_formA$delta
# rssr_rbrbp_formA <- result_rbrbp_formA$rssr
# 
# # Plot convergence
# plot(result_rbrbp_formA$all_rssr, type = "l", xlab = "Iteration", ylab = "RSSR",
#      main = "RBRBP Convergence (Form A)")

## ----eval=FALSE---------------------------------------------------------------
# # Load data
# data(example)
# y <- example$y
# bigt <- length(y)
# z <- matrix(1, nrow = bigt, ncol = 1)
# q <- 1
# m <- 3
# 
# # Form B Restrictions: delta_1 = delta_3, delta_2 = delta_4
# R <- matrix(c(1, 0, -1, 0, 0, 1, 0, -1), nrow = 2, ncol = 4, byrow = TRUE)
# rr <- c(0, 0)
# 
# # Initial estimation (Form B)
# result_initial <- est2(y, z, q = q, m = m, bigt = bigt, trm = 0.10, R = R, rr = rr)
# T_init <- result_initial$dx
# 
# cat("Initial break dates:", T_init, "\n")

## ----eval=FALSE---------------------------------------------------------------
# # Form A Restrictions: same as Form B
# S <- matrix(c(1, 0, 0, 1, 1, 0, 0, 1), nrow = 4, ncol = 2, byrow = TRUE)
# ss <- matrix(0, nrow = 4, ncol = 1)
# 
# # Initial estimation (Form A)
# result_initial_formA <- est(y, z, q = q, m = m, bigt = bigt, trm = 0.10,
#                             S = S, ss = ss)
# T_init <- result_initial_formA$dx
# 
# cat("Initial break dates:", T_init, "\n")

## ----eval=FALSE---------------------------------------------------------------
# # Run BRBP algorithm (Form B)
# result_brbp <- brbp(
#   y = y, z = z, q = q, m = m, bigt = bigt, trm = 0.10,
#   T_init = T_init,
#   B = 20,              # Number of bootstrap replications
#   R = R, rr = rr,
#   verbose = FALSE
# )
# 
# # Run RBRBP algorithm (Form B)
# result_rbrbp <- brbp_residual(
#   y = y, z = z, q = q, m = m, bigt = bigt, trm = 0.10,
#   T_init = T_init,
#   B = 20,
#   R = R, rr = rr,
#   verbose = FALSE
# )
# 
# # Compare results
# cat("Initial break dates:", T_init, "\n")
# cat("BRBP break dates:", result_brbp$dx, "\n")
# cat("RBRBP break dates:", result_rbrbp$dx, "\n")
# cat("BRBP RSSR:", result_brbp$rssr, "\n")
# cat("RBRBP RSSR:", result_rbrbp$rssr, "\n")

## ----eval=FALSE---------------------------------------------------------------
# # Run BRBP algorithm (Form A)
# result_brbp_formA <- brbp(
#   y = y, z = z, q = q, m = m, bigt = bigt, trm = 0.10,
#   T_init = T_init,
#   B = 20,
#   S = S, ss = as.vector(ss),
#   R = NULL, rr = NULL,
#   verbose = FALSE
# )
# 
# # Run RBRBP algorithm (Form A)
# result_rbrbp_formA <- brbp_residual(
#   y = y, z = z, q = q, m = m, bigt = bigt, trm = 0.10,
#   T_init = T_init,
#   B = 20,
#   S = S, ss = as.vector(ss),
#   R = NULL, rr = NULL,
#   verbose = FALSE
# )
# 
# # Compare results
# cat("Initial break dates:", T_init, "\n")
# cat("BRBP (Form A) break dates:", result_brbp_formA$dx, "\n")
# cat("RBRBP (Form A) break dates:", result_rbrbp_formA$dx, "\n")
# cat("BRBP (Form A) RSSR:", result_brbp_formA$rssr, "\n")
# cat("RBRBP (Form A) RSSR:", result_rbrbp_formA$rssr, "\n")

## -----------------------------------------------------------------------------
sessionInfo()