the main components of every time series: Trend, Seasonality, Noise and Other. (大势,小周期,噪音和其他)。
以下一段具体解释。
The trend may be linear or nonlinear, and we may want to investigate its magnitude. The
seasonality pattern may be either additive or multiplicative. In the first case, the seasonal
change has the same absolute size no matter what the magnitude of the current baseline of
the series is; in the latter case, the seasonal change has the same relative size compared
with the current magnitude of the series. Noise (i.e., some form of random variation) is
almost always part of a time series. Finding ways to reduce the noise in the data is usually
a significant part of the analysis process. Finally, “other” includes anything else that we
may observe in a time series, such as particular significant changes in overall behavior,
special outliers, missing data—anything remarkable at all.
然后就是: Description, Prediction, and Control.
Smoothing
窗口平滑,加权窗口平滑,高斯加权窗口平滑。
以上都有缺点:1, 无法评估效果,不能重复。2,由于窗口问题,不能接近真实值。3,对于范围外的点没法算,也就是不能预测。
克服上述缺点的方法:exponential smoothing or Holt–Winters method
https://gist.github.com/andrequeiroz/5888967
# Holt-Winters algorithms to forecasting # Coded in Python 2 by: Andre Queiroz # Description: This module contains three exponential smoothing algorithms. They are Holt‘s linear trend method and Holt-Winters seasonal methods (additive and multiplicative). # References: # Hyndman, R. J.; Athanasopoulos, G. (2013) Forecasting: principles and practice. http://otexts.com/fpp/. Accessed on 07/03/2013. # Byrd, R. H.; Lu, P.; Nocedal, J. A Limited Memory Algorithm for Bound Constrained Optimization, (1995), SIAM Journal on Scientific and Statistical Computing, 16, 5, pp. 1190-1208. from sys import exit from math import sqrt from numpy import array from scipy.optimize import fmin_l_bfgs_b def RMSE(params, *args): Y = args[0] type = args[1] rmse = 0 if type == ‘linear‘: alpha, beta = params a = [Y[0]] b = [Y[1] - Y[0]] y = [a[0] + b[0]] for i in range(len(Y)): a.append(alpha * Y[i] + (1 - alpha) * (a[i] + b[i])) b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i]) y.append(a[i + 1] + b[i + 1]) else: alpha, beta, gamma = params m = args[2] a = [sum(Y[0:m]) / float(m)] b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2] if type == ‘additive‘: s = [Y[i] - a[0] for i in range(m)] y = [a[0] + b[0] + s[0]] for i in range(len(Y)): a.append(alpha * (Y[i] - s[i]) + (1 - alpha) * (a[i] + b[i])) b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i]) s.append(gamma * (Y[i] - a[i] - b[i]) + (1 - gamma) * s[i]) y.append(a[i + 1] + b[i + 1] + s[i + 1]) elif type == ‘multiplicative‘: s = [Y[i] / a[0] for i in range(m)] y = [(a[0] + b[0]) * s[0]] for i in range(len(Y)): a.append(alpha * (Y[i] / s[i]) + (1 - alpha) * (a[i] + b[i])) b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i]) s.append(gamma * (Y[i] / (a[i] + b[i])) + (1 - gamma) * s[i]) y.append(a[i + 1] + b[i + 1] + s[i + 1]) else: exit(‘Type must be either linear, additive or multiplicative‘) rmse = sqrt(sum([(m - n) ** 2 for m, n in zip(Y, y[:-1])]) / len(Y)) return rmse def linear(x, fc, alpha = None, beta = None): Y = x[:] if (alpha == None or beta == None): initial_values = array([0.3, 0.1]) boundaries = [(0, 1), (0, 1)] type = ‘linear‘ parameters = fmin_l_bfgs_b(RMSE, x0 = initial_values, args = (Y, type), bounds = boundaries, approx_grad = True) alpha, beta = parameters[0] a = [Y[0]] b = [Y[1] - Y[0]] y = [a[0] + b[0]] rmse = 0 for i in range(len(Y) + fc): if i == len(Y): Y.append(a[-1] + b[-1]) a.append(alpha * Y[i] + (1 - alpha) * (a[i] + b[i])) b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i]) y.append(a[i + 1] + b[i + 1]) rmse = sqrt(sum([(m - n) ** 2 for m, n in zip(Y[:-fc], y[:-fc - 1])]) / len(Y[:-fc])) return Y[-fc:], alpha, beta, rmse def additive(x, m, fc, alpha = None, beta = None, gamma = None): Y = x[:] if (alpha == None or beta == None or gamma == None): initial_values = array([0.3, 0.1, 0.1]) boundaries = [(0, 1), (0, 1), (0, 1)] type = ‘additive‘ parameters = fmin_l_bfgs_b(RMSE, x0 = initial_values, args = (Y, type, m), bounds = boundaries, approx_grad = True) alpha, beta, gamma = parameters[0] a = [sum(Y[0:m]) / float(m)] b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2] s = [Y[i] - a[0] for i in range(m)] y = [a[0] + b[0] + s[0]] rmse = 0 for i in range(len(Y) + fc): if i == len(Y): Y.append(a[-1] + b[-1] + s[-m]) a.append(alpha * (Y[i] - s[i]) + (1 - alpha) * (a[i] + b[i])) b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i]) s.append(gamma * (Y[i] - a[i] - b[i]) + (1 - gamma) * s[i]) y.append(a[i + 1] + b[i + 1] + s[i + 1]) rmse = sqrt(sum([(m - n) ** 2 for m, n in zip(Y[:-fc], y[:-fc - 1])]) / len(Y[:-fc])) return Y[-fc:], alpha, beta, gamma, rmse def multiplicative(x, m, fc, alpha = None, beta = None, gamma = None): Y = x[:] if (alpha == None or beta == None or gamma == None): initial_values = array([0.0, 1.0, 0.0]) boundaries = [(0, 1), (0, 1), (0, 1)] type = ‘multiplicative‘ parameters = fmin_l_bfgs_b(RMSE, x0 = initial_values, args = (Y, type, m), bounds = boundaries, approx_grad = True) alpha, beta, gamma = parameters[0] a = [sum(Y[0:m]) / float(m)] b = [(sum(Y[m:2 * m]) - sum(Y[0:m])) / m ** 2] s = [Y[i] / a[0] for i in range(m)] y = [(a[0] + b[0]) * s[0]] rmse = 0 for i in range(len(Y) + fc): if i == len(Y): Y.append((a[-1] + b[-1]) * s[-m]) a.append(alpha * (Y[i] / s[i]) + (1 - alpha) * (a[i] + b[i])) b.append(beta * (a[i + 1] - a[i]) + (1 - beta) * b[i]) s.append(gamma * (Y[i] / (a[i] + b[i])) + (1 - gamma) * s[i]) y.append((a[i + 1] + b[i + 1]) * s[i + 1]) rmse = sqrt(sum([(m - n) ** 2 for m, n in zip(Y[:-fc], y[:-fc - 1])]) / len(Y[:-fc])) return Y[-fc:], alpha, beta, gamma, rmse
http://adorio-research.org/wordpress/?p=1230
def holtwinters(y, alpha, beta, gamma, c, debug=True): """ y - time series data. alpha , beta, gamma - exponential smoothing coefficients for level, trend, seasonal components. c - extrapolated future data points. 4 quarterly 7 weekly. 12 monthly The length of y must be a an integer multiple (> 2) of c. """ #Compute initial b and intercept using the first two complete c periods. ylen =len(y) if ylen % c !=0: return None fc =float(c) ybar2 =sum([y[i] for i in range(c, 2 * c)])/ fc ybar1 =sum([y[i] for i in range(c)]) / fc b0 =(ybar2 - ybar1) / fc if debug: print "b0 = ", b0 #Compute for the level estimate a0 using b0 above. tbar =sum(i for i in range(1, c+1)) / fc print tbar a0 =ybar1 - b0 * tbar if debug: print "a0 = ", a0 #Compute for initial indices I =[y[i] / (a0 + (i+1) * b0) for i in range(0, ylen)] if debug: print "Initial indices = ", I S=[0] * (ylen+ c) for i in range(c): S[i] =(I[i] + I[i+c]) / 2.0 #Normalize so S[i] for i in [0, c) will add to c. tS =c / sum([S[i] for i in range(c)]) for i in range(c): S[i] *=tS if debug: print "S[",i,"]=", S[i] # Holt - winters proper ... if debug: print "Use Holt Winters formulae" F =[0] * (ylen+ c) At =a0 Bt =b0 for i in range(ylen): Atm1 =At Btm1 =Bt At =alpha * y[i] / S[i] + (1.0-alpha) * (Atm1 + Btm1) Bt =beta * (At - Atm1) + (1- beta) * Btm1 S[i+c] =gamma * y[i] / At + (1.0 - gamma) * S[i] F[i]=(a0 + b0 * (i+1)) * S[i] print "i=", i+1, "y=", y[i], "S=", S[i], "Atm1=", Atm1, "Btm1=",Btm1, "At=", At, "Bt=", Bt, "S[i+c]=", S[i+c], "F=", F[i] print i,y[i], F[i] #Forecast for next c periods: for m in range(c): print "forecast:", (At + Bt* (m+1))* S[ylen + m] # the time-series data. y =[146, 96, 59, 133, 192, 127, 79, 186, 272, 155, 98, 219] holtwinters(y, 0.2, 0.1, 0.05, 4)