possibly correct implementation of fourier_kexp_2

tests/holt_winters.py

@@ -60,20 +60,20 @@ def modified_holt_winters(y, s_init=None, h=1, alpha=0.8, beta=0.1, gamma=0.1, o
     return preds
 
 
-m = 24
+m = 24*6
 y_hat = []
 seasonal_state = []
 s_init = data['all'].iloc[:m].copy().values
-for i in range(600):
-  preds, s = modified_holt_winters(data['all'].iloc[:1+i].copy().values, h=24, alpha=0.5, beta=0.5, gamma=0.5, m=24, return_s=True, s_init=s_init)
+for i in range(1600):
+  preds, s = modified_holt_winters(data['all'].iloc[:1+i].copy().values, h=m, alpha=0.01, beta=0.01, gamma=0.01, omega=0.1, m=m, return_s=True, s_init=s_init)
   y_hat.append(preds)
   seasonal_state.append(s)
 
 y_hat = np.vstack(y_hat)
 seasonal_state = np.vstack(seasonal_state)
 
 
-#ts_animation(data['all'].values, y_hat, 300, labels=['ground truth', 'predictions'])
+#ts_animation(data['all'].values, y_hat, 1600, labels=['ground truth', 'predictions'])
 
 
 # Fourier Holt-Winters: in this Holt-Winters forecasters the states are the Fourier coefficients
@@ -111,10 +111,10 @@ def __init__(self, h=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3):
         self.m = m
         Ps_past = []
         Ps_future = []
+        self.n_harmonics = np.minimum(n_harmonics, m // 2)
         for i in range(m):
             t = np.arange(m)+i
             t_f = np.arange(m+h)+i
-            self.n_harmonics = np.minimum(n_harmonics, m // 2)
             P = get_basis(t, m, self.n_harmonics)
             Ps_past.append(P)
             P_f = get_basis(t_f, m, self.n_harmonics)
@@ -124,6 +124,7 @@ def __init__(self, h=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3):
         self.Ps_past = Ps_past
         self.coeffs = None
         self.eps = None
+        self.last_1sa_preds = 0
     def predict(self, y, return_coeffs=False, start_from=0):
         if self.coeffs is None:
             coeffs = np.zeros(2 * self.n_harmonics + 1)
@@ -136,17 +137,20 @@ def predict(self, y, return_coeffs=False, start_from=0):
             preds = [[y[start_from]+eps]]
 
         coeffs_t_history = []
-        for i in range(1 + start_from, len(y)+1):
+        for i in range(start_from, len(y)):
             P = self.Ps_past[(i)%self.m]
-            P_f = self.Ps_future[(i)%self.m]
-            start = np.maximum(0, i-self.m)
-            w = y[start:i].copy()
-            eps = self.omega * (w[-1] - preds[-1][0]) + (1-self.omega) * eps
-            w[-1] = w[-1] -eps
+            P_f = self.Ps_future[(i-self.h)%self.m]
+            start = np.maximum(0, i-self.m+1)
+
+            w = y[start:i+1].copy()
+            #eps = self.omega * (w[-1] - preds[-1][0]) + (1-self.omega) * eps
+            eps = self.omega * (w[-1] - self.last_1sa_preds) + (1-self.omega) * eps
+            #w[-1] = w[-1] -eps
             coeffs_t = P[-len(w):,:].T@w
             coeffs = self.alpha*coeffs_t + (1-self.alpha)*coeffs
-
-            preds.append((P_f@coeffs).ravel() + eps*(self.h-np.arange(self.h))**2/self.h**2)
+            last_preds = (P_f@coeffs).ravel()
+            self.last_1sa_preds = last_preds[0]
+            preds.append((last_preds + eps*(self.h-np.arange(self.h))**2/self.h**2).ravel() )
             if return_coeffs:
                 coeffs_t_history.append(coeffs)
         self.coeffs = coeffs
@@ -156,9 +160,9 @@ def predict(self, y, return_coeffs=False, start_from=0):
         return np.vstack(preds[1:])
 
     def __getstate__(self):
-        return (self.coeffs, self.eps)
+        return (self.coeffs, self.eps,self.last_1sa_preds)
     def __setstate__(self, state):
-        self.coeffs, self.eps = state
+        self.coeffs, self.eps, self.last_1sa_preds = state
 
 from pyforecaster.forecaster import LinearForecaster
 from pyforecaster.formatter import Formatter
@@ -175,13 +179,14 @@ def __setstate__(self, state):
 ts_animation(y.iloc[n_tr:n_tr+n_te,0].values, y_hat_lin.values, 2000, labels=['ground truth', 'predictions'])
 
 
-m = 24*6
+m = int(24*6)
 
 
-y_hat, coeffs_t = fourier_hw(h=24*6, alpha=0.9, omega=0.9, n_harmonics=20, m=m).predict(data['all'].values[:2100], return_coeffs=True)
-"""
+y_hat, coeffs_t = fourier_hw(h=24*6, alpha=0.5, omega=0.1, n_harmonics=80, m=m).predict(data['all'].values[:2100], return_coeffs=True)
 ts_animation(x.iloc[:, 0].values, y_hat, 2000, labels=['ground truth', 'predictions'])
+
 ts_animation(x.iloc[:, 0].values, coeffs_t[:, :], 2000, labels=['ground truth', 'coeffs'])
+"""
 
 
 mae = lambda x, y: np.mean(np.abs(x-y))
@@ -197,6 +202,72 @@ def __setstate__(self, state):
 from filterpy.common import Q_discrete_white_noise
 
 
+class Fourier_kexp(fourier_hw):
+    def __init__(self, h=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3):
+        super().__init__(h, alpha, m, omega, n_harmonics)
+        n_harmonics = np.minimum(n_harmonics, m // 2)
+
+        coeffs = np.zeros(2 * n_harmonics + 1)
+
+        # precompute basis over all possible periods
+        my_filter = KalmanFilter(dim_x=n_harmonics * 2 + 1, dim_z=n_harmonics * 2 + 1)
+        my_filter.x = coeffs.copy()
+
+        my_filter.F = np.eye(n_harmonics * 2 + 1)
+
+        my_filter.H = np.eye(n_harmonics * 2 + 1)
+        my_filter.P = np.eye(n_harmonics * 2 + 1) * 1000
+        my_filter.R = np.eye(n_harmonics * 2 + 1) * 0.01
+        my_filter.Q = 0.01
+        self.filter = my_filter
+        self.n_harmonics = n_harmonics
+        self.coeffs_t_history = []
+    def predict(self, y, return_coeffs=False, start_from=0):
+        if self.coeffs is None:
+            coeffs = np.zeros(2 * self.n_harmonics + 1)
+            coeffs[0] = y[0]*np.sqrt(m)
+            eps = 0
+            preds = [[0]]
+        else:
+            eps = self.eps
+            coeffs = self.coeffs
+            preds = [[y[start_from]+eps]]
+
+        coeffs_t_history = []
+        for i in range(start_from, len(y)):
+            P = self.Ps_past[(i)%self.m]
+            P_f = self.Ps_future[(i-self.h)%self.m]
+            start = np.maximum(0, i-self.m+1)
+
+            w = y[start:i+1].copy()
+            #eps = self.omega * (w[-1] - preds[-1][0]) + (1-self.omega) * eps
+            eps = self.omega * (w[-1] - self.last_1sa_preds) + (1-self.omega) * eps
+            #w[-1] = w[-1] -eps
+            coeffs_t = P[-len(w):,:].T@w
+            self.filter.predict()
+            self.filter.update(coeffs_t)
+            coeffs = self.filter.x
+
+            last_preds = (P_f@coeffs).ravel()
+            self.last_1sa_preds = last_preds[0]
+            preds.append((last_preds + eps*(self.h-np.arange(self.h))**2/self.h**2).ravel() )
+            self.coeffs_t_history.append(coeffs)
+            if i % m == 0 and i >= m:
+                self.filter.Q = np.corrcoef(np.vstack(self.coeffs_t_history).T)
+                self.filter.R = np.corrcoef(np.vstack(self.coeffs_t_history).T) * 0.01
+                self.filter.P = np.eye(self.n_harmonics * 2 + 1) * 1000
+
+        self.coeffs = coeffs
+        self.eps = eps
+        if return_coeffs:
+            return np.vstack(preds[1:]), np.vstack(self.coeffs_t_history)
+        return np.vstack(preds[1:])
+
+    def __getstate__(self):
+        return (self.coeffs, self.eps,self.last_1sa_preds)
+    def __setstate__(self, state):
+        self.coeffs, self.eps, self.last_1sa_preds = state
+
 def fourier_kexp(y, h=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3, return_coeffs=False):
     """
     :param y:
@@ -236,7 +307,7 @@ def fourier_kexp(y, h=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3, return_coef
         Ps_future.append(P_f)
 
     eps = 0
-    for i in range(1, len(y)-1):
+    for i in range(1, len(y)):
         P = Ps_past[i%m]
         P_f = Ps_future[i%m]
         start = np.maximum(0, i-m)
@@ -270,9 +341,10 @@ def fourier_kexp(y, h=1, alpha=0.8, m=24, omega=0.99, n_harmonics=3, return_coef
 ts_animation(data_te, y_hat_kexp, 2000, labels=['ground truth', 'predictions'])
 ts_animation(data_te, coeffs_t, 2000, labels=['ground truth', 'predictions'])
 
+y_hat_kexp, coeffs_t = Fourier_kexp(h=24*6, alpha=0.9, omega=0.9, n_harmonics=20, m=m).predict(data_te, return_coeffs=True)
+ts_animation(data_te, y_hat_kexp, 2000, labels=['ground truth', 'predictions'])
 
-
-def fourier_kexp_2(y, h=1, m=24, omega=0.99, n_harmonics=3, return_coeffs=False):
+def fourier_kexp_2(y, h=1, n_predictors=4,  m_max=24, omega=0.9, alpha=0.9, n_harmonics=3, return_coeffs=False):
     """
     :param y:
     :param h:
@@ -281,13 +353,13 @@ def fourier_kexp_2(y, h=1, m=24, omega=0.99, n_harmonics=3, return_coeffs=False)
     :return:
     """
 
-    n_predictors = 4
-    n_harmonics = np.minimum(n_harmonics, m // 2)
 
-    coeffs = np.zeros(2*n_harmonics+1)
+    n_harmonics = np.minimum(n_harmonics, m_max // 2)
+    ms = np.linspace(1, m_max, n_predictors+1).astype(int)[1:]
+
     preds = [[0]]
     coeffs_t_history = []
-    coeffs[0] = y[0]*np.sqrt(m)
+
 
     # precompute basis over all possible periods
     my_filter = KalmanFilter(dim_x=n_predictors, dim_z=n_predictors)
@@ -299,7 +371,7 @@ def fourier_kexp_2(y, h=1, m=24, omega=0.99, n_harmonics=3, return_coeffs=False)
     my_filter.R = np.eye(n_predictors)
     my_filter.Q = 0.1
 
-    models = [fourier_hw(h=h, alpha=0.9, omega=0.9, n_harmonics=n_harmonics, m=h*(i+1)) for i in range(n_predictors)]
+    models = [Fourier_kexp(h=h, alpha=alpha, omega=omega, n_harmonics=n_harmonics, m=ms[i]) for i in range(n_predictors)]
     states = [models[j].__getstate__() for j in range(n_predictors)]
 
 
@@ -309,14 +381,11 @@ def fourier_kexp_2(y, h=1, m=24, omega=0.99, n_harmonics=3, return_coeffs=False)
         [models[j].__setstate__(states[j]) for j in range(n_predictors)]
         preds_t = [models[j].predict(y[:i+1].copy(), return_coeffs=False, start_from=i).ravel() for j in range(n_predictors)]
         for j in range(n_predictors):
-            try:
-                preds_models[j].append(preds_t[j])
-            except:
-                print(3)
+            preds_models[j].append(preds_t[j])
         if i>=h:
             for j in range(n_predictors):
                 preds_models[j].pop(0)
-        # averate last point error over different prediction times for all the models
+            # average last point error over different prediction times for all the models
             avg_err = [np.mean([np.abs(p[-(k+1)]-y[i]) for k, p in enumerate(preds_models[j])]) for j in range(n_predictors)]
             coeffs_t = np.exp(-np.array(avg_err)) / np.exp(-np.array(avg_err)).sum()
         else:
@@ -334,6 +403,7 @@ def fourier_kexp_2(y, h=1, m=24, omega=0.99, n_harmonics=3, return_coeffs=False)
         if i%m==0 and i>=m:
             my_filter.Q = np.corrcoef(np.vstack(coeffs_t_history).T)
             my_filter.R = np.corrcoef(np.vstack(coeffs_t_history).T)*10
+            my_filter.P = np.eye(n_predictors) * 10
 
     if return_coeffs:
         return np.vstack(preds[1:]), np.vstack(coeffs_t_history)
@@ -347,8 +417,8 @@ def fourier_kexp_2(y, h=1, m=24, omega=0.99, n_harmonics=3, return_coeffs=False)
 n_tr = 10000
 n_te = 2800
 
-m = 24*6*7
-y_hat_kexp, coeffs_t = fourier_kexp_2(x.iloc[n_tr:n_tr+n_te,0].values, h=24*6, omega=0.9, n_harmonics=20, m=m, return_coeffs=True)
+m = 24*6*3
+y_hat_kexp, coeffs_t = fourier_kexp_2(x.iloc[n_tr:n_tr+n_te,0].values, h=24*6, omega=0.99, alpha=0.8, n_harmonics=50, m_max=m, n_predictors=3, return_coeffs=True)
 ts_animation(x.iloc[n_tr:n_tr+n_te,0].values, y_hat_kexp, 2000, labels=['ground truth', 'predictions'])
 ts_animation(x.iloc[n_tr:n_tr+n_te,0].values, coeffs_t, 2000, labels=['ground truth', 'predictions'])