ウィンタースの季節性線形指数平滑法による予測 with Excel
Step 0イントロダクション
2009年第1四半期から2014年第4四半期の6年にわたる統計「コンビニエンスストア販売額等(合計)」のデータです(単位:百万円)。[データの出所:商業動態統計「コンビニエンスストア商品別販売額等及び前年(度,同期,同月)比」経済産業省ウェブサイト]
のちの工程でグラフによって推移をざっと見てみると,固定的なトレンド性および季節性がうかがえるところです。かかる点において,先々の見通しに関してはウィンタース法(ただし乗法)を用いることができそうです。以下,これによって将来の4四半期(2015年第1四半期~第4四半期)の販売額を確認することを目的とします。
作表・作図の上での制約
- 「ソルバー」アドインが利用可能な環境でないとおそらく困難をともないます。
作表・作図にあたり参考にした書籍・WebページおよびPDF文書
- [Web]Triple Exponential Smoothing ―"NIST/SEMATECH e-handbook Statistical Methods"(2015.10 閲覧)
- スピーロス マクリダキス・スティーブン C フィールライト(1995)『計画策定と意思決定のための予測手法入門』加藤五郎訳, 同友館, pp.64-79.
- [Web]Exponential smoothing ―"Wikipedia" (en)(2015.10 閲覧)
- [PDF]The Holt-Winters Forecasting Method, Methodology of the Monthly Index of Services Annex B: ―"Office for National Statistic"(2015.10 閲覧)
- [VideoHostingService]Forecasting in Excel using the Holt-Winters technique [id:scmprofrutgers] ―"YouTube"(2015.10 閲覧)
- [Web]Holt-Winters seasonal method ―"Otexts"(2015.10 閲覧)
「予測シート」の作成 “without” Excel2016-
Step 1初期値の計算(1)
元表がきれいな(?)うちにピボットテーブルで必要な集計・計算をおこないます。
元表のデータ範囲の任意のセルをアクティブにして「ピボットテーブルの作成」ダイアログを呼び出し,あたらしいシートにピボット表をつくります。
Step 2初期値の計算(2)
「列」に“年”の変数を,「行」に“四半期”の変数を割り当て,「値」に“実測値(ここでは販売額)”を割り当てます。
“実測値”が「合計」以外になった場合は(「データの個数」など),これを「合計」に修正します。
同じくピボット表・行見出し側の「総計」を変更します。これをアクティブにした上でフィールドの設定ボタンから
集計方法を「平均」に変更します。
Step 3初期値の計算(3)
ピボット表のままでは以降の諸計算をおこないづらいので,ここでは下図赤色囲みの領域を表のすぐ下にコピーして以降の処理を進めるものとします。
あたらしくコピーした方の表の最下行(年平均)の最も古い2サイクル(ここでは2年)について,それらの差をとり(順序注意),期数(ここでは4つの四半期=4)で除します。これを期間当初の1期あたりの平均的なトレンドとみなし,のち工程で利用していきます(ここではこれを便宜上,b0と呼びます)。
セルH16 | =(C16-B16)/4 |
---|
表のさらに下方で,あらたに下図のような計算をおこないます。実測値を年平均で除し,四半期ごとの趨勢を指数化します(対象期間分すべてのセルにコピー)。
セルB18 | =B12/B$16 |
---|
この表の右方で上で求めた値から四半期別に平均をとり,必要な分だけコピーします(ここでは残り3行)。これを対象期間における各四半期の(販売額の)平均的な指数とみなし,こちらものちの工程で利用します(ここではこれを便宜上,c0と呼びます)。
セルH18 | =AVERAGE(B18:G18) |
---|
Step 4予測シートの準備
ピボット表のシートを離れ,当初のシートに戻ります。
シートの各所に,下図のような見出しを用意します。なおここで表中に登場する列見出しには,次の意味を与えています(見出しなので任意のもので可)。
- [alpha,beta,gamma] パラメータ(ウエイト)
- [s] 基準値(季節調整済平滑実測値)
- [b] トレンドの平滑値
- [c] 季節値(季節指数)の平滑値
- [F] 予測値(当期を境に先後)
- [RMSE] 平均二乗誤差
Step 5初期値の入力
Step 3で導いたb0(下図)を[b]の初期値とします。
ピボット表から当該セルをコピーし
元シートの下図の位置(2サイクル目[ここでは2年目]が始まる直前)に値のみ貼り付けます。
同様に,c0(下図)を[c]の初期値とします。
ピボット表から当該セルをコピーし
元シートの側の下図の範囲(最初の4四半期分)に値のみを貼り付けます。
残る[s]の初期値には,最初の4四半期分の平均値を充てるものとします(2サイクル目[ここでは2年目]が始まる直前に)。
セルE8 | =AVERAGE(D5:D8) |
---|
Step 6初期パラメータの入力
alpha,beta,gammaの初期値を設定します。
これらは後段でソルバーにてfitする値を別途探索していくので,この時点では0<alpha<1,0<beta<1,0<gamma<1の条件を満たす文字通りテキトーな任意の値を充てておきます(ここではすべて0.1とします)。
Step 7基準値の計算
以下,各列の計算を処理します。最初に[s]列を計算します。
[s]は下式によって求めます。
L: 1サイクルの期数(この例では4)。以下に同じ。
つまり,当期の季節調整済み実測値にalpha,トレンドの平滑値を加算した前期の基準値(季節調整値)に1-alphaの係数(ウエイト)をそれぞれ掛けて合成します。
これを2サイクル目の最初(t=5)のセルに入力すると,下表・下図のようになります。
セルE9 | =$E$2*D9/G5+(1-$E$2)*(E8+F8) |
---|
Step 8トレンドの計算
次に[b]列を計算します。
[b]は下式によって求めます。
つまり当期・前期の基準値の差にbeta,前期のトレンドの平滑値に1-betaの係数(ウエイト)をそれぞれ掛けて合成します。
これを2サイクル目の最初(t=5)のセルに入力すると,下表・下図のようになります。
セルF9 | =$F$2*(E9-E8)+(1-$F$2)*F8 |
---|
Step 9季節指数の計算
つづいて[c]列を計算します。
[c]は下式によって求めます。
つまり当期実測値・基準値の比にgamma,前サイクルの同じ期(この例では「四半期」)期のトレンドの平滑値に1-gammaの係数(ウエイト)をそれぞれ掛けて合成します。
これを2サイクル目の最初(t=5)のセルに入力すると,下表・下図のようになります。
セルG9 | =$G$2*D9/E9+(1-$G$2)*G5 |
---|
Step 10主要3式のコピー
Step 7-9の3式を表の最下行までコピーします。
Step 11予測式へ過去の変数の値をあてはめる
[s], [b], [c]の各値をもとに各期における予測値を求めます。ここでいう予測値はグラフ描画の観点から2種に分けるものとします(実測値の揃った前段"←F"において導く「予測式を通した値」と,それ以降の予測値"F→")。
[F]は下式によって求めます。
m: 予測したいのはベースとする[s], [b]の値のm期あと。
つまり,ウィンタース法での予測値は,当期の基準値にm期分のトレンド平滑値を加算した値を季節指数で乗じて求められます。
まずは「←F」列を埋めることとします。この列ではすべて最も小さなmの値(m=1)で予測値が計算できるので,その点をふまえて式を作成すると,下表・下図のようになります。またこの式を表の最下行までコピーします。
セルH9 | =(E8+F8*1)*G5 |
---|
Step 12誤差の計算
予測値が実測値とできるだけfitするようなパラメータが欲しいので,その判断のよりどころとしての誤差を求めます。ここではRMSEと呼ばれる指標に依存することとします(これが小さな値を返すほどあてはまりがよいと判断)。
この値を求めるための手続きのひとつとして,「(x-F)^2」列で実測値と予測値との差の2乗値を求めます(上式Σよりあとの部分を計算)。
セルJ9 | =(D9-H9)^2 |
---|
セルJ2 | =SQRT(AVERAGE(J9:J28)) |
---|
ふたつ目として,見出し「RMSE」直下に下表・下図のような式を入力します(残りの部分を計算)※。
※ より古い期のfitの具合はあまり重視しないなど,筋論から言えば本来は“どの程度の範囲を評価の対象とするか”によってAverage関数の引数(データ範囲)も検討される必要もありそうです。しかしここでは機械的に,予測値を出せる全期間を対象にするものとして進めます。
Step 13次サイクルの予測値の計算(1)
本来の意味での予測値(いまだ実測値の存在しない期間の予測値)を導くにあたり,t,Y,Qといった項目を1サイクル分延長します。
実測値と予測値をグラフ上でスムースにつなげるため,t,Y,Qを除く他の列の最下行と同じ「F→」列のセルに,当期の実測値を転記します。
Step 14次サイクルの予測値の計算(2)
再びStep 11の式を使って,「F→」列の予測値を求めます。
具体的には,下表・下図の通りです。
セルI29 | =($E$28+$F$28*(A29-$A$28))*G25 |
---|
またこの計算式を先に延長したセルの分だけ下方にコピーしておきます。
Step 15グラフの作成
Q※, x, ←F, F→の4列を選択し,折れ線グラフを作成します。
※ 作例では内容が四半期番号のみなので,実際には別途作業列などにラベル用の列を作成した方が良さげです。この場合,2013より前のバージョンではConcatenateで,2013以降ではフラッシュフィルでY-Q列の文字列を連結させると簡単です。
上記のシートの状態からグラフを作成すると,下図のようになります。このグラフにいう青線が実測値,オレンジの線が予測式をとおした各期の値(いわば,過去の時点での予測値),グレーの線が予測値となります。
したがって,青い線とオレンジの線が(RMSEの計算において指定した範囲で)重なるほどfitしたモデルと見ることができます。
Step 16パラメータの探索[ソルバー](1)
よりfitするパラメータ(があるかどうかを含めて)を探すため,alpha,beta,gammaの値をいろいろと変更してみます。アタリをつけながら手作業でこれをおこなうのは多くの場合難儀するので,ここではソルバー機能を利用します。
「ソルバーのパラメータ」ダイアログを呼び出し,「目的セルの設定」にRMSEを求めたセルを,「目標値」に最小値を,「変数セルの変更」にalpha,beta,gamma一連のセルを指定します(RMSEが最小になるalpha,beta,gammaの組み合わせを探したい)。
「制約条件の対象」の追加ボタンをクリックします(alpha,beta,gammaは0より大きくて1より小さな値でなければならない)。
「制約条件の追加」ダイアログで3つのパラメータに関する制約を2個つくります。
- alpha,beta,gammaは0より大きい
- セル参照:alpha ~ gamma 値のセル範囲
- 大小関係:>=
- 制約条件:0
なお当該条件から厳密にいえば大小関係は“>”ですがソルバーの仕様もありここでは原則を無視します(次に同じ)。
- alpha,beta,gammaは1より小さい
- セル参照:alpha ~ gamma 値のセル範囲
- 大小関係:<=
- 制約条件:1
すべての設定をおえて,制約条件が下図のようになっているか確認します。
Step 17パラメータの探索[ソルバー](2)
ダイアログの解決ボタンをクリックして解析を実行します。
(しばらくののち)「ソルバーの結果」ダイアログが表示されたらOKを返します(シート上のパラメータが変更されています。これによってシートも再計算されています)。
Step 18Tips1―初期値による結果の違い
解析前と解析後のグラフを並べてみると,モデルが実測値とよりfitしたかどうかがよくわかります。
なお,fitのしかたは初期値によっても異なります。下のグラフはこの例示での初期値で作成したもので,
さらに下のグラフは別の初期値を用意して作成したものです。この場合,本来の意味での予測値も当然変わります。
Step 19Tips2―予測期間の延長
1サイクル(作例では4四半期)を超えて予測値を出したい場合,計算式が直近のサイクルの[c]値を巡回して参照できるように 元式に若干の調整を挟みます。
以下,その一例です。
I29 | =IF(MOD(A29-$A$28, 4)<>0, ($E$28+$F$28*(A29-$A$28))*OFFSET($G$24, MOD(A29-$A$28, 4), 0), ($E$28+$F$28*(A29-$A$28))*OFFSET($G$24, 4, 0)) |
---|
Step 20その後の実測値
最後に参考として,この例示で用いた実測値以降は2015年第2四半期まで下図のように販売額が公表されています(2015年10月時点。青字の箇所)。
その他の参照
このサイトの関連How-toです。
メインサイト「ひとりマーケティングのためのデータ分析」の予測ツールに関するHow-toです。
同じくトレンド分析ツールに関するHow-toです。