diff --git a/tests/holt_winters.py b/tests/holt_winters.py index 8291ca2..815e7b9 100644 --- a/tests/holt_winters.py +++ b/tests/holt_winters.py @@ -60,12 +60,12 @@ 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) @@ -73,7 +73,7 @@ def modified_holt_winters(y, s_init=None, h=1, alpha=0.8, beta=0.1, gamma=0.1, o 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'])