# Psilocybin Research Database - Bai-Perron structural-break analysis
# Source: https://psilocybin-research.com/msca-impact/
# Complete years through 2025; visible, non-false-positive catalogue records.
# Snapshot SHA-256: 756f970ec11384e2b82ac95f6f7a042c56b3eab80cd8b042a7edf2f0299f1415
# Includes the annual counts: no account, database download or credentials needed.
# Run this entire file in RStudio, source("psilocybin-bai-perron-structural-breaks.R"), or Rscript psilocybin-bai-perron-structural-breaks.R.
# The default snapshot reproduces this download. To use current database counts,
# change "snapshot" to "live" below (internet required). Results may then change.

data_source <- Sys.getenv("PSILOCYBIN_DATA_SOURCE", "snapshot")
snapshot_data <- read.csv(text = paste(c(
  "year,count", "1958,5", "1959,11", "1960,11", "1961,14", "1962,15",
  "1963,20", "1964,10", "1965,10", "1966,9", "1967,12", "1968,15",
  "1969,13", "1970,12", "1971,9", "1972,3", "1973,15", "1974,15",
  "1975,9", "1976,7", "1977,9", "1978,10", "1979,7", "1980,9",
  "1981,16", "1982,8", "1983,4", "1984,7", "1985,6", "1986,8",
  "1987,6", "1988,6", "1989,6", "1990,7", "1991,1", "1992,4",
  "1993,6", "1994,6", "1995,1", "1996,5", "1997,4", "1998,15",
  "1999,18", "2000,17", "2001,8", "2002,17", "2003,18", "2004,12",
  "2005,20", "2006,25", "2007,16", "2008,24", "2009,24", "2010,16",
  "2011,22", "2012,24", "2013,37", "2014,39", "2015,41", "2016,65",
  "2017,73", "2018,106", "2019,114", "2020,170", "2021,304", "2022,423",
  "2023,550", "2024,653", "2025,885"
), collapse = "\n"), colClasses = c("integer", "integer"))
last_complete_year <- 2025L
if (identical(data_source, "snapshot")) {
  annual_counts <- snapshot_data
} else if (identical(data_source, "live")) {
  csv_file <- tempfile(fileext = ".csv")
  tryCatch({
    download.file("https://psilocybin-research.com/msca-impact/analysis-code.php?format=csv", csv_file, method = "libcurl", quiet = TRUE,
                  headers = c("User-Agent" = "Mozilla/5.0 (compatible; PsilocybinResearch-R/1.0)"))
    annual_counts <- read.csv(csv_file, colClasses = c("integer", "integer"))
  }, error = function(e) {
    stop("Could not download annual counts. Use data_source = 'snapshot' to run offline.\n", conditionMessage(e))
  }, finally = unlink(csv_file))
  last_complete_year <- as.integer(format(Sys.time(), "%Y", tz = "UTC")) - 1L
} else {
  stop("data_source must be 'snapshot' or 'live'.")
}
stopifnot(all(c("year", "count") %in% names(annual_counts)))
annual_counts <- annual_counts[complete.cases(annual_counts) &
                               annual_counts$year <= last_complete_year & annual_counts$count > 0L, ]
annual_counts <- annual_counts[order(annual_counts$year), ]
if (anyDuplicated(annual_counts$year)) stop("Annual counts must have unique years.")
cat("Source: Psilocybin Research Database\n")
cat("Data:", data_source, "| complete years through", last_complete_year, "\n")
if (identical(data_source, "snapshot")) cat("Snapshot SHA-256: 756f970ec11384e2b82ac95f6f7a042c56b3eab80cd8b042a7edf2f0299f1415\n")

# Descriptive break dating: global least-squares segmentation with BIC.
# This base-R implementation reproduces the website's explicit BIC convention.
# BIC = n*log(RSS/n) + (2*segments + breaks)*log(n).
# Minimum segment: max(5, ceiling(0.15*n)); search up to five breaks.
# Break years mark the START of a new segment. These are estimates, not exact dates.
# This script does not assign p-values or confidence intervals to the dates.
# Reference: Bai & Perron (2003), https://doi.org/10.1002/jae.659

estimate_breaks <- function(data, max_breaks = 5L, trim = 0.15) {
  n <- nrow(data)
  if (n < 5L) stop("At least five positive-count annual observations are required.")
  h <- max(5L, ceiling(trim * n))
  max_segments <- min(max_breaks + 1L, n %/% h)
  x <- data$year - min(data$year)
  y <- log(data$count)
  rss <- matrix(Inf, n, n)
  for (start in seq_len(n - h + 1L)) {
    for (end in seq.int(start + h - 1L, n)) {
      ix <- seq.int(start, end)
      fit <- lm.fit(cbind(1, x[ix]), y[ix])
      rss[start, end] <- sum(fit$residuals^2)
    }
  }
  floor_rss <- max(1e-12, rss[1L, n] * 1e-12)
  cost <- matrix(Inf, max_segments + 1L, n + 1L)
  parent <- matrix(NA_integer_, max_segments + 1L, n + 1L)
  cost[1L, 1L] <- 0
  candidates <- data.frame(segments = seq_len(max_segments), rss = NA_real_, bic = NA_real_)
  for (segments in seq_len(max_segments)) {
    for (end in seq.int(segments * h, n)) {
      starts <- seq.int((segments - 1L) * h, end - h)
      costs <- cost[segments, starts + 1L] + rss[cbind(starts + 1L, rep(end, length(starts)))]
      best <- which.min(costs)
      if (is.finite(costs[best])) {
        cost[segments + 1L, end + 1L] <- costs[best]
        parent[segments + 1L, end + 1L] <- starts[best]
      }
    }
    candidates$rss[segments] <- cost[segments + 1L, n + 1L]
    candidates$bic[segments] <- n * log(max(candidates$rss[segments], floor_rss) / n) +
      (2 * segments + segments - 1L) * log(n)
  }
  best <- 1L
  for (i in seq_len(nrow(candidates))) {
    if (candidates$bic[i] < candidates$bic[best] - 1e-9) best <- i
  }
  cuts <- integer(0)
  end <- n
  for (segments in seq.int(best, 1L)) {
    start <- parent[segments + 1L, end + 1L]
    if (start > 0L) cuts <- c(start, cuts)
    end <- start
  }
  list(break_years = data$year[cuts + 1L], observations = n,
       minimum_segment = h, max_breaks = max_breaks,
       bic = candidates$bic[best], delta_bic = candidates$bic[1L] - candidates$bic[best],
       candidates = candidates)
}
analysis_result <- estimate_breaks(annual_counts)
print(analysis_result)
if (interactive()) {
  plot(annual_counts$year, log(annual_counts$count), type = "b", pch = 19,
       col = "#476f91", xlab = "Publication year", ylab = "log(catalogue records)",
       main = "BIC-selected structural breaks")
  abline(v = analysis_result$break_years, col = "#c8943f", lty = 2, lwd = 2)
}

breakpoint_result <- analysis_result

# Formal structural-change tests on the same annual-count snapshot.
# Formal Bai-Perron tests for structural change, not separate p-values per date.
# Validated against mbreaks 1.0.1. First run installs that version if necessary.
# The package is pure R; its ggplot2 dependency is installed through CRAN if absent.
# Internet is needed for initial package setup, not for the included count snapshot.
# Reference: Bai & Perron (2003), https://doi.org/10.1002/jae.659

ensure_mbreaks <- function(version = "1.0.1") {
  installed <- tryCatch(as.character(packageVersion("mbreaks")), error = function(e) "")
  if (identical(installed, version)) return(invisible(TRUE))
  if (isNamespaceLoaded("mbreaks")) {
    stop("A different mbreaks version is already loaded. Start a fresh R session and run this file again.")
  }
  # Keep the pinned analysis package separate from other R projects.
  library_dir <- Sys.getenv("PSILOCYBIN_R_LIBRARY",
                            file.path(path.expand("~"), "R", "psilocybin-research", paste0("R-", getRversion())))
  dir.create(library_dir, recursive = TRUE, showWarnings = FALSE)
  if (file.access(library_dir, 2L) != 0L) stop("Cannot create a writable R package folder in your user directory.")
  .libPaths(c(library_dir, .libPaths()))
  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    install.packages("ggplot2", repos = "https://cloud.r-project.org", lib = library_dir)
  }
  package_file <- tempfile(fileext = ".tar.gz")
  on.exit(unlink(package_file), add = TRUE)
  candidates <- c(
    sprintf("https://cloud.r-project.org/src/contrib/mbreaks_%s.tar.gz", version),
    sprintf("https://cloud.r-project.org/src/contrib/Archive/mbreaks/mbreaks_%s.tar.gz", version)
  )
  downloaded <- FALSE
  for (address in candidates) {
    status <- tryCatch(suppressWarnings(download.file(address, package_file, quiet = TRUE, mode = "wb")),
                       error = function(e) 1L)
    if (identical(status, 0L)) { downloaded <- TRUE; break }
  }
  if (!downloaded) stop("Could not download mbreaks 1.0.1. Check your internet connection and run this file again.")
  install.packages(package_file, repos = NULL, type = "source", lib = library_dir)
  if (!identical(as.character(packageVersion("mbreaks")), version)) {
    stop("mbreaks installation did not complete. Review the installation messages above.")
  }
}
ensure_mbreaks()

# More conservative testing settings than descriptive BIC dating:
# 20% trimming rounded DOWN, up to 3 breaks; changing intercept and time slope.
# Segment-specific HAC covariance, quadratic-spectral kernel, Andrews AR(1)
# bandwidth and VAR(1) prewhitening. Assumes short-range serial dependence.
# Finite annual samples and catalogue coverage require cautious interpretation.
d <- annual_counts
if (nrow(d) < 40L || any(diff(d$year) != 1L)) {
  stop("Formal tests require at least 40 consecutive positive-count complete years.")
}
d$log_count <- log(d$count)
d$time <- (d$year - min(d$year)) / (max(d$year) - min(d$year))
global_tests <- mbreaks::dotest("log_count", "time", data = d,
                              m = 3L, eps1 = 0.20, robust = 1L, prewhit = 1L,
                              hetdat = 1L, hetvar = 1L)
sequential_tests <- mbreaks::doseqtests("log_count", "time", data = d,
                                      m = 3L, eps1 = 0.20, robust = 1L, prewhit = 1L,
                                      hetdat = 1L, hetvar = 1L)
p_bound <- function(statistic, critical) {
  alpha <- as.numeric(rownames(critical)) / 100
  for (i in order(alpha)) {
    if (statistic > critical[i, 1L]) return(paste("p <", format(alpha[i], scientific = FALSE)))
  }
  "p >= 0.1"
}
udmax <- max(as.numeric(global_tests$ftest))
analysis_result <- list(
  observations = nrow(d), first_year = min(d$year), last_year = max(d$year),
  minimum_segment = floor(0.20 * nrow(d)), max_breaks = 3L,
  udmax = udmax, udmax_p_bound = p_bound(udmax, global_tests$cv_Dmax),
  supf = as.numeric(global_tests$ftest),
  supf_p_bounds = vapply(seq_len(3L), function(k) p_bound(global_tests$ftest[k, 1L], global_tests$cv_supF[, k, drop = FALSE]), character(1L)),
  sequential = as.numeric(sequential_tests$supfl),
  sequential_p_bounds = vapply(seq_len(3L), function(k) {
    value <- sequential_tests$supfl[k, 1L]
    if (value <= 0) NA_character_ else p_bound(value, sequential_tests$cv[, k, drop = FALSE])
  }, character(1L)),
  mbreaks_version = as.character(packageVersion("mbreaks"))
)
print(analysis_result)
print(global_tests)
print(sequential_tests)
cat("\nP-values are bounds from published asymptotic critical tables, not exact probabilities.\n")
cat("The overall test does not give a separate p-value or confidence interval for each break year.\n")
cat("UDmax uses the package's standard double-max critical table; searching only three alternatives is conservative.\n")
cat("A zero sequential statistic in mbreaks can mean no feasible extra split; the website shows that as unavailable.\n")

formal_test_result <- analysis_result
analysis_result <- breakpoint_result