Chapter Four, Time As a Variable: Time-Series Analysis

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
View Code

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)
View Code
原文地址:https://www.cnblogs.com/hluo/p/4065148.html