Metadata-Version: 2.4
Name: online-svr
Version: 1.0.0
Summary: Online SVR for time series forecasting with budgeted working set and adaptive tolerance.
Author-email: Ademir Neto <ademir.bsn@gmail.com>
License: MIT
Project-URL: Repository, https://github.com/ademirNeto/online-svr
Keywords: svr,time-series,forecasting,online-learning,budgeted-learning
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy
Requires-Dist: scikit-learn
Requires-Dist: scipy
Requires-Dist: statsmodels
Provides-Extra: dev
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: pytest-cov; extra == "dev"
Requires-Dist: jupyter; extra == "dev"
Requires-Dist: nbconvert; extra == "dev"
Dynamic: license-file

# online-svr

A Python library for time series forecasting with **online Support Vector Regression**, featuring a scikit-learn compatible API and support for adaptive tolerance triggering with a budgeted working set.

## Installation

```bash
pip install online-svr
```

## Overview

`online-svr` implements the algorithm described in "[article submitted to *Evolving Systems*]":

- **Budgeted working set**: The number of support vectors is kept within a budget `M` (default 1500) by removing the vector with the smallest RKHS-norm contribution `α_i² K(x_i, x_i)`.
- **Adaptive tolerance**: The update trigger uses a threshold `κ · s_t`, where `s_t` is an EWMA of the absolute residuals (RiskMetrics, λ=0.94).
- **Representation selection**: Automatically chooses the best representation (levels vs. first differences) for your series.
- **One-step-ahead forecasting**: Prequential interface for online evaluation and learning.

## Quick Start

```python
import numpy as np
from online_svr import OnlineSVRForecaster

rng = np.random.default_rng(42)
y = np.cumsum(rng.standard_normal(5000))

f = OnlineSVRForecaster(M=1500, kappa=2.0, lam=0.94)
f.fit(y[:4000])

preds = []
for t in range(4000, 4500):
    tail = y[t - f.p_ - 1: t]
    yhat = f.predict_update(y[t], tail)
    preds.append(yhat)

print(f"Predicted {len(preds)} steps.")
```

## Key Classes

| Class | Description |
|---|---|
| `BudgetedOnlineSVR` | Online SVR with budgeted working set and adaptive tolerance. Core algorithm. |
| `OnlineSVRForecaster` | Time-series wrapper: lag selection (PACF), representation auto-selection, one-step-ahead interface. |

## Metrics

Evaluation functions for forecasting:

| Function | Description |
|---|---|
| `theil_u2` | Theil's U2 (comparison to random walk; U2 < 1 is better). |
| `mase` | Mean Absolute Scaled Error (scale-invariant). |
| `mape` | Mean Absolute Percentage Error (use with caution on series crossing zero). |
| `pocid` | Prediction of Change In Direction (fraction of correct direction predictions). |

```python
from online_svr import theil_u2, mase, pocid

y_true = y[4000:4500]
y_naive = np.concatenate([[y[3999]], y[3999:4499]])

print(f"U2 vs random walk: {theil_u2(y_true, np.array(preds)):.3f}")
print(f"POCID: {pocid(y_true, np.array(preds)):.3f}")
```

## How It Works

### BudgetedOnlineSVR

1. **`fit_initial(X, y, resid_scale_init)`** — Offline phase. Fits an SVR on the training data, initializes the working set with support vectors, and prunes to the budget.
2. **`predict_update(x, y_true)`** — Prequential step: predict on `x`, then decide whether to add `x` to the working set based on the residual vs. `κ · s_t`. If added, refit the SVR and prune if necessary.

### OnlineSVRForecaster

1. **`fit(y_history, val_fraction=1/8)`** — Offline phase. Selects significant lags via PACF, chooses the representation (levels or first differences) by validation MSE, and seeds the online model.
2. **`predict_update(y_new, history_tail)`** — One-step-ahead prediction and update in one call. Returns the forecast made *before* seeing `y_new`.

### Representation Selection

- **`representation="levels"`**: Model the series directly.
- **`representation="naive-residual"`**: Model the first differences and add them to the last observation.
- **`representation="auto"`** (default): Fit both and choose the one with lower validation MSE.

The "naive-residual" mode is particularly useful for non-stationary series and automatically competes with a random-walk baseline.

## Architecture & Flow

### Component Responsibilities

| Component | Role |
|---|---|
| `BudgetedOnlineSVR` | Core algorithm: SVR fitting + working-set budget + adaptive trigger |
| `OnlineSVRForecaster` | Time-series wrapper: lag selection, representation choice, prequential API |
| `StandardScaler` | Feature normalization (fitted on training data, applied online) |
| `SVR` (sklearn) | Underlying estimator (refit on the working set at each trigger) |

### Prequential Loop (Sequence Diagram)

```
┌──────────────────────────────────────────────────────────────────┐
│ Offline Phase: fit(y_train)                                       │
│   1. Select lags via PACF                                         │
│   2. Choose representation (levels vs naive-residual)             │
│   3. Grid search on validation block (temporal holdout)           │
│   4. Seed working set = support vectors from training             │
│   5. Initialize s_t = median(|residuals|) on validation           │
└──────────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌──────────────────────────────────────────────────────────────────┐
│ Online Phase: for t in y_test:                                    │
│   1.  predict(x_t)         →  ŷ_t                                │
│   2.  e_t = |y_t - ŷ_t|                                            │
│   3.  if e_t > κ·s_t:                                                │
│         add (x_t, y_t) to working set                              │
│         refit SVR on working set                                   │
│         if len(working set) > M:                                   │
│           prune: keep top-M by α_i² K(x_i, x_i)                  │
│   4.  s_t ← √(λ·s_t² + (1-λ)·e_t²)   [EWMA, RiskMetrics]         │
└──────────────────────────────────────────────────────────────────┘
```

**Key invariants:**
- Working set is **always bounded** by `M` (no memory blowup).
- The model is **never** refit on the full history, only on the working set.
- `s_t` is the only state that persists across steps (along with the working set).

## Use Cases

The three examples below use synthetic data so you can run them without external files.

### Use Case 1: Stationary Series (Easy)

A stationary series has a stable mean and no trend. The SVR can learn the noise pattern and beat a random-walk baseline.

```python
import numpy as np
from online_svr import OnlineSVRForecaster, theil_u2, mape

rng = np.random.default_rng(42)
# AR(1) around zero
y = np.zeros(10000)
for i in range(1, len(y)):
    y[i] = 0.7 * y[i-1] + rng.standard_normal()

f = OnlineSVRForecaster(M=500, kappa=2.0, representation="levels", 
                        max_lags=5, init_window=7000)
f.fit(y[:8500])

preds = [f.predict_update(y[t], y[t - f.p_ - 1: t]) for t in range(8500, 10000)]
y_true = y[8500:10000]
print(f"U2 = {theil_u2(y_true, np.array(preds)):.3f}")  # expect ~0.5-0.8
print(f"MAPE = {mape(y_true, np.array(preds)):.2%}")
```

**Lesson:** With a stationary target, the SVR learns the autocorrelation and beats the random walk.

### Use Case 2: Non-Stationary Series (Intermediate)

Random walks (or trending data) are hard to forecast. The `naive-residual` representation models the increment `y_t - y_{t-1}`, which is often stationary even when `y_t` is not.

```python
import numpy as np
from online_svr import OnlineSVRForecaster, theil_u2

rng = np.random.default_rng(7)
y = np.cumsum(rng.standard_normal(10000))  # random walk

f = OnlineSVRForecaster(M=1000, kappa=2.0, representation="auto", 
                        max_lags=5, init_window=7000)
f.fit(y[:8500])

preds = [f.predict_update(y[t], y[t - f.p_ - 1: t]) for t in range(8500, 10000)]
y_true = y[8500:10000]
print(f"Representation chosen: {f.rep_}")
print(f"U2 = {theil_u2(y_true, np.array(preds)):.3f}")  # ~1.0 (close to naive)
```

**Lesson:** `representation="auto"` will pick `naive-residual` for non-stationary data. U2 near 1.0 means "as good as the random walk" — which is often the ceiling for genuine random walks.

### Use Case 3: Seasonal Series with Outliers (Advanced)

Real-world series (electricity load, monthly sales) have seasonality and occasional spikes. Two strategies help:
- Increase `kappa` to avoid reacting to outliers as new support vectors.
- Increase `lam` to make `s_t` more inertial, smoothing single-event shocks.

```python
import numpy as np
from online_svr import OnlineSVRForecaster, theil_u2

rng = np.random.default_rng(99)
n = 10000
t = np.arange(n)
# Trend + weekly seasonality + rare outliers
y = 0.01 * t + 2 * np.sin(2 * np.pi * t / 50) + 0.3 * rng.standard_normal(n)
# Inject 2% outliers
outlier_idx = rng.choice(n, size=int(0.02 * n), replace=False)
y[outlier_idx] += rng.choice([-1, 1], size=len(outlier_idx)) * 15

f = OnlineSVRForecaster(
    M=1200, kappa=3.0, lam=0.96,  # more tolerant + smoother EWMA
    representation="auto", max_lags=8, init_window=7000,
)
f.fit(y[:8500])

preds = [f.predict_update(y[t_], y[t_ - f.p_ - 1: t_]) for t_ in range(8500, n)]
y_true = y[8500:n]
print(f"U2 = {theil_u2(y_true, np.array(preds)):.3f}")  # expect <1.0
print(f"Refits: {f.model_.n_refits}")  # show trigger count
```

**Lesson:** Outliers inflate residuals and may trigger spurious updates. A larger `kappa` filters them out; a larger `lam` makes `s_t` absorb shocks more slowly.

## Parameter Tuning Guide

### `OnlineSVRForecaster`

| Parameter | Default | Range | Effect | When to increase |
|---|---|---|---|---|
| `M` | 1500 | 500–5000 | Max size of the working set | More memory available; series with complex patterns |
| `kappa` | 2.0 | 1.0–5.0 | Multiplier of the adaptive tolerance | Noisy series (reduces spurious refits); clean series (decrease for more updates) |
| `lam` | 0.94 | 0.90–0.98 | EWMA smoothing (RiskMetrics) | 0.90 = more responsive to recent residuals; 0.98 = more inertial |
| `representation` | `"auto"` | `"levels"`, `"naive-residual"`, `"auto"` | Target for the model | Use `"auto"` if unsure; use `"levels"` only when you know the series is stationary |
| `max_lags` | 12 | 3–24 | Max lags evaluated via PACF | Slow dynamics (increase); fast dynamics (decrease) |
| `init_window` | 8000 | N/2–N | Max training samples used to seed the working set | Long series (irrelevant cap); short series (decrease) |

### Practical Recommendations

- **Start with defaults.** Observe `f.model_.n_refits` (how often the model updates) and U2.
- **If U2 >> 1.0:** Increase `kappa` (reduces false triggers). The series may also be inherently unpredictable — check MASE.
- **If U2 << 1.0 but you're hitting the budget:** Increase `M` (more memory for diverse patterns) or reduce `max_lags` (fewer features).
- **If the working set is always full:** Consider lowering `kappa` (more selective refits) or check for data quality issues (duplicates, outliers).
- **For long production runs:** Monitor `len(f.model_.Wx)` and `f.model_.n_refits`. A stable working set size is a sign that the model has converged on the relevant patterns.

### When NOT to use this library

- **Series shorter than ~8000 observations:** The default `init_window` won't fit. Use `init_window = N // 2` manually or pick a different model.
- **Tabular regression:** SVR is not the right tool; use `sklearn.ensemble` or gradient boosting.
- **High-frequency data with strong short-term autocorrelation:** Try `representation="levels"` with `max_lags>=20`.

## Limitations

- Series must have at least `init_window + max_lags + 5` observations for `fit` to succeed.
- Only PSD kernels (linear, rbf) are supported.
- For very short series or extreme outliers, the adaptive tolerance may need manual tuning.
- The grid search in `fit` evaluates 160 hyperparameter combinations by default; on long series this can take a few seconds.

## References

If you use this library in your research, please cite the corresponding article:

> Neto, A. (202X). Online support vector regression with budgeted working set and adaptive tolerance. *Evolving Systems*, ...

## License

MIT
