///|
/// A regularly ordered numeric time series.
pub struct CausalTimeSeries {
  time : Array[Int]
  values : Array[Double]
  frequency : Int
  name : String
}

///|
/// Time-series diagnostics used before an interrupted-series analysis.
pub struct TimeSeriesProfile {
  observations : Int
  missing : Int
  mean : Double
  standard_deviation : Double
  first_difference_mean : Double
  autocorrelation_lag1 : Double
  trend_slope : Double
  seasonal_period : Int
  passes : Bool
}

///|
/// Forecast result with residual scale and horizon labels.
pub struct ForecastResult {
  time : Array[Int]
  forecast : Array[Double]
  lower : Array[Double]
  upper : Array[Double]
  residual_standard_error : Double
  model_order : Int
}

///|
/// Interrupted time-series estimate.
pub struct InterruptedEffect {
  level_change : Double
  slope_change : Double
  standard_error : Double
  pre_slope : Double
  post_slope : Double
  pre_observations : Int
  post_observations : Int
  passes : Bool
}

///|
/// Candidate change point diagnostic.
pub struct ChangePoint {
  index : Int
  statistic : Double
  left_mean : Double
  right_mean : Double
  detected : Bool
}

///|
/// Newey-West long-run variance summary.
pub struct LongRunVariance {
  variance : Double
  standard_error : Double
  bandwidth : Int
  autocovariances : Array[Double]
}

///|
fn ts_mean(values : Array[Double]) -> Double {
  mean_or(values, 0.0)
}

///|
fn ts_time_mean(time : Array[Int]) -> Double {
  if time.length() == 0 {
    0.0
  } else {
    let mut total = 0.0
    for value in time {
      total += value.to_double()
    }
    total / time.length().to_double()
  }
}

///|
/// Creates a validated time series.
pub fn causal_time_series(
  time : Array[Int],
  values : Array[Double],
  frequency? : Int = 1,
  name? : String = "outcome",
) -> CausalTimeSeries {
  let n = time.length().min(values.length())
  {
    time: time[:n].to_owned(),
    values: values[:n].to_owned(),
    frequency: if frequency > 0 {
      frequency
    } else {
      1
    },
    name,
  }
}

///|
/// Returns a profile for a time series.
pub fn time_series_profile(
  series : CausalTimeSeries,
  seasonal_period? : Int = 1,
) -> TimeSeriesProfile {
  let missing = missing_count(series.values)
  let observed = series.values.filter(fn(value) { is_finite(value) })
  let differences = first_difference(series.values)
  let slope = linear_trend_slope(series.time, series.values)
  let autocorrelation = autocorrelation_at(series.values, 1)
  {
    observations: series.values.length(),
    missing,
    mean: ts_mean(observed),
    standard_deviation: std_dev(observed),
    first_difference_mean: ts_mean(differences),
    autocorrelation_lag1: autocorrelation,
    trend_slope: slope,
    seasonal_period: if seasonal_period > 0 {
      seasonal_period
    } else {
      1
    },
    passes: missing == 0 && series.values.length() >= 4,
  }
}

///|
/// Returns a lagged vector aligned with the original index.
pub fn lag_series(
  values : Array[Double],
  lag : Int,
  fill? : Double = 0.0,
) -> Array[Double] {
  let result = Array::make(values.length(), fill)
  if lag <= 0 {
    return values.copy()
  }
  for i in lag.. Array[Double] {
  let result = Array::make(values.length(), fill)
  if lead <= 0 {
    return values.copy()
  }
  if lead < values.length() {
    for i in 0..<(values.length() - lead) {
      result[i] = values[i + lead]
    }
  }
  result
}

///|
/// Computes first differences with a configurable initial value.
pub fn first_difference(
  values : Array[Double],
  initial? : Double = 0.0,
) -> Array[Double] {
  let result = Array::new(capacity=values.length())
  if values.length() == 0 {
    return result
  }
  result.push(values[0] - initial)
  for i in 1.. Array[Double] {
  if lag <= 0 {
    return values.copy()
  }
  let result = Array::new(capacity=values.length())
  for i in lag.. Array[Double] {
  let result = Array::new(capacity=values.length())
  if values.length() == 0 {
    return result
  }
  let width = if window > 0 { window } else { 1 }
  for i in 0.. width { i + 1 - width } else { 0 }
    let end = i + 1
    let slice = values[start:end].to_owned()
    result.push(ts_mean(slice))
  }
  result
}

///|
/// Computes a trailing moving sum.
pub fn moving_sum(values : Array[Double], window : Int) -> Array[Double] {
  let result = Array::new(capacity=values.length())
  let width = if window > 0 { window } else { 1 }
  let mut running = 0.0
  for i in 0..= width {
      running -= values[i - width]
    }
    result.push(running)
  }
  result
}

///|
/// Computes exponentially weighted moving averages.
pub fn exponential_moving_average(
  values : Array[Double],
  alpha : Double,
) -> Array[Double] {
  let result = Array::new(capacity=values.length())
  if values.length() == 0 {
    return result
  }
  let weight = clamp(alpha, 1.0e-6, 1.0)
  let mut level = values[0]
  result.push(level)
  for value in values[1:] {
    level = weight * value + (1.0 - weight) * level
    result.push(level)
  }
  result
}

///|
/// Computes rolling standard deviations.
pub fn rolling_standard_deviation(
  values : Array[Double],
  window : Int,
) -> Array[Double] {
  let result = Array::new(capacity=values.length())
  let width = if window > 0 { window } else { 1 }
  for i in 0.. width { i + 1 - width } else { 0 }
    let end = i + 1
    result.push(std_dev(values[start:end].to_owned()))
  }
  result
}

///|
/// Computes the lag-k sample autocorrelation.
pub fn autocorrelation_at(values : Array[Double], lag : Int) -> Double {
  if lag <= 0 || values.length() <= lag {
    return 0.0
  }
  let average = mean(values)
  let mut numerator = 0.0
  let mut denominator = 0.0
  for value in values {
    let difference = value - average
    denominator += difference * difference
  }
  for i in lag.. Array[Double] {
  let result = Array::new(
    capacity=if maximum_lag > 0 { maximum_lag } else { 0 },
  )
  for lag in 1..<=maximum_lag {
    result.push(autocorrelation_at(values, lag))
  }
  result
}

///|
/// Computes a least-squares linear trend slope.
pub fn linear_trend_slope(time : Array[Int], values : Array[Double]) -> Double {
  let n = time.length().min(values.length())
  if n < 2 {
    return 0.0
  }
  let x_mean = ts_time_mean(time[:n].to_owned())
  let y_mean = ts_mean(values[:n].to_owned())
  let mut numerator = 0.0
  let mut denominator = 0.0
  for i in 0.. Double {
  mean_or(values, 0.0) - linear_trend_slope(time, values) * ts_time_mean(time)
}

///|
/// Fits an autoregressive model of order one by ordinary least squares.
pub fn fit_ar1(values : Array[Double]) -> Array[Double] {
  if values.length() < 2 {
    return [0.0, ts_mean(values)]
  }
  let lagged = values[:values.length() - 1].to_owned()
  let current = values[1:].to_owned()
  let slope = covariance(lagged, current) / variance(lagged)
  let intercept = mean(current) - slope * mean(lagged)
  [intercept, if is_finite(slope) { slope } else { 0.0 }]
}

///|
/// Fits an autoregressive model with several lags using a regularized normal equation.
pub fn fit_ar(
  values : Array[Double],
  order : Int,
  ridge? : Double = 1.0e-6,
) -> Array[Double] {
  let width = if order > 0 { order } else { 1 }
  if values.length() <= width {
    return Array::make(width + 1, 0.0)
  }
  let rows : Array[Array[Double]] = Array::new(capacity=values.length() - width)
  let target = Array::new(capacity=values.length() - width)
  for i in width.. ForecastResult {
  let steps = if horizon > 0 { horizon } else { 0 }
  let history = values.copy()
  let forecasts : Array[Double] = Array::new(capacity=steps)
  let order = if coefficients.length() > 1 {
    coefficients.length() - 1
  } else {
    1
  }
  for _ in 0.. 0 {
      coefficients[0]
    } else {
      0.0
    }
    for lag in 1..<=order {
      let index = history.length() - lag
      if index >= 0 && lag < coefficients.length() {
        prediction += coefficients[lag] * history[index]
      }
    }
    forecasts.push(prediction)
    history.push(prediction)
  }
  let residual_scale = if values.length() > order {
    std_dev(k_difference(values, 1))
  } else {
    0.0
  }
  let time = Array::new(capacity=steps)
  let lower = Array::new(capacity=steps)
  let upper = Array::new(capacity=steps)
  for i in 0.. Array[Double] {
  let result = values.copy()
  let n = result.length()
  for i in 0..= 0 && !is_finite(result[left]) {
        left -= 1
      }
      let mut right = i + 1
      while right < n && !is_finite(result[right]) {
        right += 1
      }
      if left >= 0 && right < n {
        result[i] = result[left] +
          (result[right] - result[left]) *
          (i - left).to_double() /
          (right - left).to_double()
      } else if left >= 0 {
        result[i] = result[left]
      } else if right < n {
        result[i] = result[right]
      } else {
        result[i] = 0.0
      }
    }
  }
  result
}

///|
/// Computes seasonal means by position within a cycle.
pub fn seasonal_indices(values : Array[Double], period : Int) -> Array[Double] {
  let width = if period > 0 { period } else { 1 }
  let result = Array::make(width, 0.0)
  let counts = Array::make(width, 0)
  for i in 0.. Array[Double] {
  let result = Array::new(capacity=values.length())
  if indices.length() == 0 {
    return values.copy()
  }
  for i in 0.. Array[Double] {
  let result = Array::new(capacity=values.length())
  if indices.length() == 0 {
    return values.copy()
  }
  for i in 0.. InterruptedEffect {
  let pre_time = Array::new()
  let pre_values = Array::new()
  let post_time = Array::new()
  let post_values = Array::new()
  for i in 0..= 3 && post_values.length() >= 3,
  }
}

///|
/// Scans a series for the largest mean-shift statistic.
pub fn detect_change_point(
  values : Array[Double],
  minimum_segment? : Int = 3,
) -> ChangePoint {
  let n = values.length()
  let minimum = if minimum_segment > 0 { minimum_segment } else { 3 }
  if n < 2 * minimum {
    return {
      index: 0,
      statistic: 0.0,
      left_mean: 0.0,
      right_mean: 0.0,
      detected: false,
    }
  }
  let mut best_index = minimum
  let mut best_statistic = 0.0
  let mut best_left = 0.0
  let mut best_right = 0.0
  for index in minimum..<(n - minimum + 1) {
    let left = values[:index].to_owned()
    let right = values[index:].to_owned()
    let difference = (mean(left) - mean(right)).abs()
    let pooled = (variance(left) / index.to_double() +
    variance(right) / (n - index).to_double()).sqrt()
    let statistic = if pooled == 0.0 { difference } else { difference / pooled }
    if statistic > best_statistic {
      best_statistic = statistic
      best_index = index
      best_left = mean(left)
      best_right = mean(right)
    }
  }
  {
    index: best_index,
    statistic: best_statistic,
    left_mean: best_left,
    right_mean: best_right,
    detected: best_statistic > 2.0,
  }
}

///|
/// Computes autocovariances through a Newey-West bandwidth.
pub fn newey_west_variance(
  values : Array[Double],
  bandwidth : Int,
) -> LongRunVariance {
  let n = values.length()
  if n < 2 {
    return {
      variance: 0.0,
      standard_error: 0.0,
      bandwidth: 0,
      autocovariances: [],
    }
  }
  let width = if bandwidth > 0 { bandwidth.min(n - 1) } else { 0 }
  let center = mean(values)
  let autocovariances : Array[Double] = Array::new(capacity=width + 1)
  for lag in 0..<=width {
    let mut covariance_value = 0.0
    for i in lag.. Int {
  if sample_size <= 2 {
    0
  } else {
    (4.0 * @math.pow(sample_size.to_double() / 100.0, 0.25))
    .round()
    .to_int()
    .min(sample_size - 1)
  }
}

///|
/// Computes a block bootstrap index sequence for dependent observations.
pub fn block_bootstrap_indices(
  sample_size : Int,
  block_size : Int,
  seed : UInt64,
) -> Array[Int] {
  let n = if sample_size > 0 { sample_size } else { 0 }
  let width = if block_size > 0 { block_size.min(n.max(1)) } else { 1 }
  let result = Array::new(capacity=n)
  let rng = RandomState::new(seed)
  while result.length() < n {
    let start = if n == 0 {
      0
    } else {
      (rng.uniform() * n.to_double()).to_int()
    }
    for offset in 0.. Array[Double] {
  let means = moving_average(values, window)
  let deviations = rolling_standard_deviation(values, window)
  let result = Array::new(capacity=values.length())
  for i in 0.. Array[Bool] {
  let scores = rolling_z_scores(values, window)
  let result = Array::new(capacity=values.length())
  for score in scores {
    result.push(!is_finite(score) || score.abs() > threshold)
  }
  result
}

///|
/// Computes a seasonal naive forecast.
pub fn seasonal_naive_forecast(
  values : Array[Double],
  period : Int,
  horizon : Int,
  start_time? : Int = 0,
) -> ForecastResult {
  let width = if period > 0 { period } else { 1 }
  let steps = if horizon > 0 { horizon } else { 0 }
  let forecast = Array::new(capacity=steps)
  for i in 0..= 0 && index < values.length() {
        values[index]
      } else {
        ts_mean(values)
      },
    )
  }
  let error = if values.length() > width {
    std_dev(k_difference(values, width))
  } else {
    0.0
  }
  let time = Array::new(capacity=steps)
  let lower = Array::new(capacity=steps)
  let upper = Array::new(capacity=steps)
  for i in 0.. Array[Double] {
  let profile = time_series_profile(series, seasonal_period=series.frequency)
  [
    profile.observations.to_double(),
    profile.missing.to_double(),
    profile.mean,
    profile.standard_deviation,
    profile.first_difference_mean,
    profile.autocorrelation_lag1,
    profile.trend_slope,
    profile.seasonal_period.to_double(),
    if profile.passes {
      1.0
    } else {
      0.0
    },
  ]
}