時間序列分解
大量時間序列的觀測樣本表現出趨勢性、季節性和隨機性,或者三者中的其一或其二。於是,我們認為每個時間序列,都可以分為三個部分的疊加
其中,T是趨勢項,S是季節項,R是隨機項。
上述公式表現了趨勢項和季節項是累加的,實際應用場景中,趨勢項和季節項可能是累乘的,時間序列可以分解為如下公式
實際應用中,隨機項R的期望為0,沒有規律,並且絕對值不大。所以在應用場景中我們往往省略掉R,R稱作噪聲。預測公式如下
或
一次指數平滑法
線性回歸算法中,每個經驗點的權重是一致的,即很早以前的經驗數據也可能對預測數據有較大的影響。很多實際場景中,未來一段時間的趨勢可能和在最近一段時間的趨勢關系更加緊密。比如小明去年數學考試成績一直不及格,今年連續多次考試90多分,預測小明下一次數學考試的成績,情理上90多分的可能性更高。采用傳統的線性回歸算法,預測結果可能是70多分。
指數平滑法認為越老的經驗數據對趨勢的影響越小。我們假定時間t的觀測值為y(t),時間t的預測值為S(t),則時間t+1的預測值S(t+1)為
a的取值范圍(0, 1),a越大,最近時間點的觀測值對預測值的影響越大。
假設我們有t個經驗數據,根據上述一次指數平滑公式,預測值S(t + n) = S(t + 1),預測值不具備趨勢。
二次指數平滑
我們對一次指數平滑值再進行指數平滑,可以獲得趨勢。二次指數平滑法的預測模型為:
式中:分別為時間t和時間t - 1的二次指數平滑值。
三次指數平滑
二次指數模型是線性的,對於非線性趨勢預測我們可以使用三次指數平滑法。公式如下
Holt-Winters算法
對於具有周期性的趨勢預測,我們可以使用Holt-Winters算法。累乘性Holt-Winters公式如下
其中,alpha,beta,gamma取值范圍為(0, 1),分別表示全局因子,趨勢因子,周期性因子中最近時間點數據對預測數據的影響程度。y為經驗數據,L為周期。
表示使用t時間點的估計值預測t+m時間點的值。
注:預測公式中I(t – L + m)應該為I(t – L + 1 + (m – 1) mod L)
計算步驟
alpha,beta,gamma,y,L,m已知。
(1)初始化S0
S0 = y0
(2)初始化b0
(3)初始化I1, I2, …, IL
(3)計算所有S,b,I
(4)根據公式預測未來值。其中,t取經驗數據最后一個時間點,t+m為預測時間點。
累加性Holt-Winters公式
Holt-Winters理解
指數平滑法與Holt-Winters不是建立在理論基礎上的,而是一種經驗法則。文章開發我們討論了時間序列的分解,Holt-Winters公式正是把時間序列分解為趨勢項和周期項。其中趨勢項為線性函數s + mb,周期項為c。這里面趨勢項與周期項考慮了指數平滑,即給不同時間點的趨勢或者周期性賦予了不同的權重。不過,這里的趨勢僅僅是線性趨勢,在帶有季節性的非線性趨勢預測中,效果可能不那么好。
Holt-Winters的Java實現
下面的代碼實在google上搜索的,預測公式邏輯有問題,沒有對季節項進行mod運算,會導致數組越界,有空再修改代碼。
package coshaho.learn; public class HoltWinters { public static double[] forecast(int[] y, double alpha, double beta, double gamma, int period, int m, boolean debug) { if (y == null) { return null; } int seasons = y.length / period; double a0 = calculateInitialLevel(y, period); double b0 = calculateInitialTrend(y, period); double[] initialSeasonalIndices = calculateSeasonalIndices(y, period, seasons); if (debug) { System.out.println(String.format( "Total observations: %d, Seasons %d, Periods %d", y.length, seasons, period)); System.out.println("Initial level value a0: " + a0); System.out.println("Initial trend value b0: " + b0); printArray("Seasonal Indices: ", initialSeasonalIndices); } double[] forecast = calculateHoltWinters(y, a0, b0, alpha, beta, gamma, initialSeasonalIndices, period, m, debug); if (debug) { printArray("Forecast", forecast); } return forecast; } private static double[] calculateHoltWinters(int[] y, double a0, double b0, double alpha, double beta, double gamma, double[] initialSeasonalIndices, int period, int m, boolean debug) { double[] St = new double[y.length]; double[] Bt = new double[y.length]; double[] It = new double[y.length]; double[] Ft = new double[y.length + m]; //Initialize base values St[1] = a0; Bt[1] = b0; for (int i = 0; i < period; i++) { It[i] = initialSeasonalIndices[i]; } Ft[m] = (St[0] + (m * Bt[0])) * It[0];//This is actually 0 since Bt[0] = 0 Ft[m + 1] = (St[1] + (m * Bt[1])) * It[1];//Forecast starts from period + 2 //Start calculations for (int i = 2; i < y.length; i++) { //Calculate overall smoothing if((i - period) >= 0) { St[i] = alpha * y[i] / It[i - period] + (1.0 - alpha) * (St[i - 1] + Bt[i - 1]); } else { St[i] = alpha * y[i] + (1.0 - alpha) * (St[i - 1] + Bt[i - 1]); } //Calculate trend smoothing Bt[i] = gamma * (St[i] - St[i - 1]) + (1 - gamma) * Bt[i - 1]; //Calculate seasonal smoothing if((i - period) >= 0) { It[i] = beta * y[i] / St[i] + (1.0 - beta) * It[i - period]; } //Calculate forecast if( ((i + m) >= period) ){ Ft[i + m] = (St[i] + (m * Bt[i])) * It[i - period + m]; } if(debug){ System.out.println(String.format( "i = %d, y = %d, S = %f, Bt = %f, It = %f, F = %f", i, y[i], St[i], Bt[i], It[i], Ft[i])); } } return Ft; } /** * See: http://robjhyndman.com/researchtips/hw-initialization/ * 1st period's average can be taken. But y[0] works better. * * @return - Initial Level value i.e. St[1] */ private static double calculateInitialLevel(int[] y, int period) { /** double sum = 0; for (int i = 0; i < period; i++) { sum += y[i]; } return sum / period; **/ return y[0]; } /** * See: http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc435.htm * * @return - Initial trend - Bt[1] */ private static double calculateInitialTrend(int[] y, int period){ double sum = 0; for (int i = 0; i < period; i++) { sum += (y[period + i] - y[i]); } return sum / (period * period); } /** * See: http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc435.htm * * @return - Seasonal Indices. */ private static double[] calculateSeasonalIndices(int[] y, int period, int seasons){ double[] seasonalAverage = new double[seasons]; double[] seasonalIndices = new double[period]; double[] averagedObservations = new double[y.length]; for (int i = 0; i < seasons; i++) { for (int j = 0; j < period; j++) { seasonalAverage[i] += y[(i * period) + j]; } seasonalAverage[i] /= period; } for (int i = 0; i < seasons; i++) { for (int j = 0; j < period; j++) { averagedObservations[(i * period) + j] = y[(i * period) + j] / seasonalAverage[i]; } } for (int i = 0; i < period; i++) { for (int j = 0; j < seasons; j++) { seasonalIndices[i] += averagedObservations[(j * period) + i]; } seasonalIndices[i] /= seasons; } return seasonalIndices; } private static void printArray(String description, double[] data){ System.out.println(String.format("******************* %s *********************", description)); for (int i = 0; i < data.length; i++) { System.out.println(data[i]); } System.out.println(String.format("*****************************************************************", description)); } }