anofox-forecast 0.15.9

Time series forecasting library
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166

#include "anofox-time/models/theta.hpp"
#include "anofox-time/models/theta_utils.hpp"
#include "anofox-time/utils/logging.hpp"
#include <cmath>
#include <numeric>
#include <algorithm>
#include <stdexcept>

namespace anofoxtime::models {

namespace {
    constexpr double kEpsilon = 1e-10;
}

Theta::Theta(int seasonal_period, double theta_param)
    : seasonal_period_(seasonal_period), theta_(theta_param),
      alpha_(0.5), level_(0.0) {
    if (seasonal_period_ < 1) {
        throw std::invalid_argument("Seasonal period must be >= 1");
    }
    if (theta_ <= 0.0) {
        throw std::invalid_argument("Theta parameter must be positive");
    }
}

std::vector<double> Theta::deseasonalize(const std::vector<double>& data) {
    std::vector<double> deseasonalized;
    theta_utils::deseasonalize(data, seasonal_period_, deseasonalized, seasonal_indices_);
    return deseasonalized;
}

std::vector<double> Theta::reseasonalize(const std::vector<double>& forecast) const {
    std::vector<double> reseasonalized;
    theta_utils::reseasonalize(forecast, seasonal_indices_, seasonal_period_, 
                               history_.size(), reseasonalized);
    return reseasonalized;
}

void Theta::fitRaw(const std::vector<double>& data) {
    if (data.empty()) {
        throw std::invalid_argument("Cannot fit Theta on empty data");
    }
    
    history_ = data;
    deseasonalized_ = deseasonalize(history_);
    
    // Use Pegels state-space formulation
    // Initial smoothed value is the first observation
    double initial_smoothed = deseasonalized_[0];
    
    // For STM (Standard Theta Method): use fixed alpha and theta
    // Alpha will be set externally or defaults to 0.5
    // Theta is set in constructor (default 2.0)
    
    states_.resize(deseasonalized_.size());
    std::vector<double> e(deseasonalized_.size());
    std::vector<double> amse(3);  // 3-step-ahead MSE
    
    // Calculate with Pegels formulation
    theta_pegels::calc(deseasonalized_, states_, 
                      theta_pegels::ModelType::STM,
                      initial_smoothed, alpha_, theta_,
                      e, amse, 3);
    
    // Store final level
    level_ = states_.back()[0];
    
    // Compute fitted values and residuals
    computeFittedValues();
    
    is_fitted_ = true;
    
    ANOFOX_INFO("Theta fitted with alpha={:.4f}, theta={:.2f}", alpha_, theta_);
}

void Theta::fit(const core::TimeSeries& ts) {
    if (ts.dimensions() != 1) {
        throw std::invalid_argument("Theta currently supports univariate series only");
    }
    
    fitRaw(ts.getValues());
}

void Theta::computeFittedValues() {
    const size_t n = deseasonalized_.size();
    fitted_.resize(n);
    residuals_.resize(n);
    
    // Fitted values from mu component
    for (size_t i = 0; i < n; ++i) {
        fitted_[i] = states_[i][4];  // mu is the forecast
        residuals_[i] = deseasonalized_[i] - fitted_[i];
    }
    
    // Reseasonalize if needed
    if (seasonal_period_ > 1 && !seasonal_indices_.empty()) {
        for (size_t i = 0; i < n; ++i) {
            size_t season_idx = i % static_cast<size_t>(seasonal_period_);
            fitted_[i] *= seasonal_indices_[season_idx];
            residuals_[i] = history_[i] - fitted_[i];
        }
    }
}

core::Forecast Theta::predict(int horizon) {
    if (!is_fitted_) {
        throw std::runtime_error("Theta::predict called before fit");
    }
    
    if (horizon <= 0) {
        return {};
    }
    
    // Generate forecast using Pegels formulation
    std::vector<double> forecast(horizon);
    theta_pegels::StateMatrix workspace(states_.size() + horizon);
    theta_pegels::forecast(states_, states_.size(), 
                          theta_pegels::ModelType::STM,
                          forecast, alpha_, theta_, workspace);
    
    // Reseasonalize
    forecast = reseasonalize(forecast);
    
    core::Forecast result;
    result.primary() = forecast;
    return result;
}

core::Forecast Theta::predictWithConfidence(int horizon, double confidence) {
    if (confidence <= 0.0 || confidence >= 1.0) {
        throw std::invalid_argument("Confidence level must be between 0 and 1");
    }
    
    auto forecast = predict(horizon);
    
    if (residuals_.empty()) {
        return forecast;
    }
    
    // Compute residual variance
    double sum_sq = 0.0;
    for (double r : residuals_) {
        sum_sq += r * r;
    }
    const double sigma = std::sqrt(sum_sq / static_cast<double>(residuals_.size()));
    
    // Normal quantile for confidence interval
    const double z = 1.96;  // Approximate 95% CI
    
    // Compute confidence intervals
    auto& lower = forecast.lowerSeries();
    auto& upper = forecast.upperSeries();
    lower.resize(horizon);
    upper.resize(horizon);
    
    for (int h = 0; h < horizon; ++h) {
        const double std_h = sigma * std::sqrt(static_cast<double>(h + 1));
        lower[h] = forecast.primary()[h] - z * std_h;
        upper[h] = forecast.primary()[h] + z * std_h;
    }
    
    return forecast;
}

} // namespace anofoxtime::models