1import numpy as np
2
3class SimpleARIMA:
4    """
5    Simplified ARIMA(p, d, q) implementation.
6    Uses ordinary least squares for parameter estimation.
7    """
8
9    def __init__(self, p=1, d=1, q=0):
10        self.p = p
11        self.d = d
12        self.q = q
13        self.ar_params = None
14        self.constant = None
15        self.original_tail = None
16
17    def _difference(self, series, d):
18        """Apply d-th order differencing."""
19        result = series.copy()
20        for _ in range(d):
21            result = np.diff(result)
22        return result
23
24    def _undifference(self, forecast, last_values):
25        """Reverse the differencing operation."""
26        result = forecast.copy()
27        for i, last_val in enumerate(reversed(last_values)):
28            cumulative = np.cumsum(result) + last_val
29            result = cumulative
30        return result
31
32    def fit(self, series):
33        """Fit ARIMA model using OLS for AR parameters."""
34        self.original_tail = series[-self.d:] if self.d > 0 else np.array([])
35
36        # Apply differencing
37        diffed = self._difference(series, self.d)
38
39        # Create lagged features for AR(p)
40        n = len(diffed)
41        if n <= self.p:
42            self.ar_params = np.zeros(self.p)
43            self.constant = np.mean(diffed)
44            return self
45
46        X = np.column_stack([
47            diffed[self.p - i - 1:n - i - 1] for i in range(self.p)
48        ])
49        y = diffed[self.p:]
50
51        # Add constant column
52        X_with_const = np.column_stack([np.ones(len(y)), X])
53
54        # OLS: params = (X^T X)^-1 X^T y
55        try:
56            params = np.linalg.lstsq(X_with_const, y, rcond=None)[0]
57            self.constant = params[0]
58            self.ar_params = params[1:]
59        except np.linalg.LinAlgError:
60            self.constant = np.mean(diffed)
61            self.ar_params = np.zeros(self.p)
62
63        self.fitted_diffed = diffed
64        return self
65
66    def predict(self, steps=10):
67        """Generate multi-step forecasts."""
68        diffed = self.fitted_diffed.copy()
69        history = list(diffed)
70        forecasts = []
71
72        for _ in range(steps):
73            # AR prediction
74            recent = np.array(history[-self.p:])[::-1]
75            pred = self.constant + np.dot(self.ar_params, recent)
76            forecasts.append(pred)
77            history.append(pred)
78
79        forecasts = np.array(forecasts)
80
81        # Undifference
82        if self.d > 0:
83            last_vals = list(self.original_tail)
84            forecasts = self._undifference(forecasts, last_vals)
85
86        return forecasts
87
88
89# Demo
90np.random.seed(42)
91t = np.arange(200)
92series = 50 + 0.1 * t + np.cumsum(np.random.randn(200) * 0.5)
93
94model = SimpleARIMA(p=2, d=1, q=0)
95model.fit(series)
96forecast = model.predict(steps=10)
97
98print(f"Last 5 actual values: {np.round(series[-5:], 2)}")
99print(f"10-step forecast: {np.round(forecast, 2)}")

1import numpy as np
2
3class ExponentialSmoothing:
4    """Simple, Holt's, and Holt-Winters exponential smoothing."""
5
6    def __init__(self, alpha=0.3, beta=None, gamma=None, seasonal_period=None):
7        self.alpha = alpha
8        self.beta = beta       # Trend smoothing
9        self.gamma = gamma     # Seasonal smoothing
10        self.seasonal_period = seasonal_period
11
12    def fit_predict(self, series, forecast_horizon=10):
13        """Fit the model and forecast."""
14        n = len(series)
15        m = self.seasonal_period
16
17        if self.beta is None and self.gamma is None:
18            # Simple Exponential Smoothing
19            return self._simple_es(series, forecast_horizon)
20        elif self.gamma is None:
21            # Holt's (trend, no seasonality)
22            return self._holts(series, forecast_horizon)
23        else:
24            # Holt-Winters (trend + seasonality)
25            return self._holt_winters(series, forecast_horizon)
26
27    def _simple_es(self, series, h):
28        level = series[0]
29        for t in range(1, len(series)):
30            level = self.alpha * series[t] + (1 - self.alpha) * level
31        return np.full(h, level)
32
33    def _holts(self, series, h):
34        level = series[0]
35        trend = series[1] - series[0]
36        for t in range(1, len(series)):
37            new_level = self.alpha * series[t] + (1 - self.alpha) * (level + trend)
38            trend = self.beta * (new_level - level) + (1 - self.beta) * trend
39            level = new_level
40        return np.array([level + (i + 1) * trend for i in range(h)])
41
42    def _holt_winters(self, series, h):
43        m = self.seasonal_period
44        # Initialize seasonal from first period
45        seasonal = np.zeros(m)
46        for i in range(m):
47            seasonal[i] = series[i] - np.mean(series[:m])
48
49        level = np.mean(series[:m])
50        trend = (np.mean(series[m:2*m]) - np.mean(series[:m])) / m if len(series) >= 2*m else 0
51
52        for t in range(m, len(series)):
53            new_level = self.alpha * (series[t] - seasonal[t % m]) + (1 - self.alpha) * (level + trend)
54            trend = self.beta * (new_level - level) + (1 - self.beta) * trend
55            seasonal[t % m] = self.gamma * (series[t] - new_level) + (1 - self.gamma) * seasonal[t % m]
56            level = new_level
57
58        forecasts = []
59        for i in range(h):
60            forecasts.append(level + (i + 1) * trend + seasonal[(len(series) + i) % m])
61        return np.array(forecasts)
62
63
64# Demo: Holt-Winters on seasonal data
65np.random.seed(42)
66t = np.arange(120)
67series = 50 + 0.2 * t + 10 * np.sin(2 * np.pi * t / 12) + np.random.randn(120) * 2
68
69model = ExponentialSmoothing(alpha=0.3, beta=0.1, gamma=0.2, seasonal_period=12)
70forecast = model.fit_predict(series, forecast_horizon=12)
71
72print(f"Last 6 values: {np.round(series[-6:], 2)}")
73print(f"12-step forecast: {np.round(forecast, 2)}")

1import numpy as np
2
3def forecast_metrics(actual, predicted):
4    """Compute common time series forecasting metrics."""
5    actual = np.array(actual)
6    predicted = np.array(predicted)
7    errors = actual - predicted
8
9    mae = np.mean(np.abs(errors))
10    rmse = np.sqrt(np.mean(errors ** 2))
11
12    # MAPE (handle zeros)
13    mask = actual != 0
14    mape = np.mean(np.abs(errors[mask] / actual[mask])) * 100 if mask.any() else np.inf
15
16    return {"MAE": mae, "RMSE": rmse, "MAPE": mape}
17
18
19def naive_forecast(train, test_len):
20    """Naive: repeat last value."""
21    return np.full(test_len, train[-1])
22
23
24def seasonal_naive(train, test_len, period=7):
25    """Seasonal naive: repeat last season."""
26    last_season = train[-period:]
27    reps = test_len // period + 1
28    return np.tile(last_season, reps)[:test_len]
29
30
31def walk_forward_validation(series, model_fn, initial_train_size=100,
32                             forecast_horizon=1, step_size=1):
33    """
34    Walk-forward (expanding window) backtesting.
35
36    model_fn: callable that takes (train_data) and returns forecast of length forecast_horizon
37    """
38    n = len(series)
39    actuals = []
40    predictions = []
41
42    for t in range(initial_train_size, n - forecast_horizon + 1, step_size):
43        train = series[:t]
44        actual = series[t:t + forecast_horizon]
45
46        pred = model_fn(train)[:forecast_horizon]
47        actuals.extend(actual)
48        predictions.extend(pred)
49
50    return forecast_metrics(actuals, predictions)
51
52
53# Demo
54np.random.seed(42)
55n = 200
56t_idx = np.arange(n)
57series = 50 + 0.1 * t_idx + 5 * np.sin(2 * np.pi * t_idx / 7) + np.random.randn(n) * 2
58
59train, test = series[:160], series[160:]
60
61# Compare baselines
62naive_pred = naive_forecast(train, len(test))
63seasonal_pred = seasonal_naive(train, len(test), period=7)
64
65print("Naive forecast:", forecast_metrics(test, naive_pred))
66print("Seasonal naive:", forecast_metrics(test, seasonal_pred))
67
68# Walk-forward validation
69simple_model = lambda train: naive_forecast(train, 1)
70wf_result = walk_forward_validation(series, simple_model, initial_train_size=100)
71print("Walk-forward (naive):", wf_result)