7.4 役に立つ予測変数

時系列データに回帰を使う場合にしばしば出てくる、役に立つ予測変数がいくつかあります。

トレンド

時系列データにトレンドがあることはよくあることです。線形トレンドは、単に\(x_{1,t}=t\)を予測変数として使うだけでモデル化できます。 \[ y_{t}= \beta_0+\beta_1t+\varepsilon_t \] ただし、\(t=1,\dots,T\)です。トレンド変数はTSLM()関数の中で特殊関数trend()を使うことで指定できます。非線形トレンドでのモデル化は7.7節で議論します。

ダミー変数

ここまでは、各予測変数の値は数値と想定してきました。しかし、予測変数が(例えば、“yes”か”no”の)2値しか取らないカテゴリー変数だったら、どうしましょう? 例えば、日次売上高を予測する際、その日が祝日か否かを考慮したい場合に、そうした変数が出てきます。祝日なら”yes”を、否なら”no”を値とする予測変数です。

この状況は、“yes”に対しては1、“no”に対しては0の値を取る「ダミー変数」を作り出すことで、重回帰モデルのフレームワーク内で扱えます。ダミー変数は「指標変数」とも言います。

データ中の外れ値を考慮するためにも、ダミー変数は使えます。外れ値を省く替わりに、ダミー変数がその効果を除去します。このケースでは、ダミー変数は外れ値なら1、そうでなければ0の値を取ります。一つの例は、特別なイベントがあったケースです。例えば、ブラジルへの到着旅行者数を予測する場合、2016年リオデジャネイロ夏季五輪の影響を考慮する必要があるでしょう。

カテゴリーが3つ以上あるカテゴリー変数の場合、複数(カテゴリー総数マイナス1個)のダミー変数を使って表わせます。ファクター型の値を持つ変数を予測変数に指定すると、TSLM()はこのケースを自動的に取り扱います。わざわざ対応するダミー変数を作り出す必要は通常ありません。

季節ダミー変数

日次データを予測していて、曜日を予測変数として考慮したいとしましょう。その場合、以下のダミー変数を作り出せます。

\(d_{1,t}\) \(d_{2,t}\) \(d_{3,t}\) \(d_{4,t}\) \(d_{5,t}\) \(d_{6,t}\)
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
0 0 0 0 0 0
1 0 0 0 0 0

7つのカテゴリーを表すのに、ダミー変数は6つだけで済むことに留意してください。7つ目のカテゴリー(このケースでは、日曜日)は、全てのダミー変数がゼロとして指定され、その影響は切片に取り込まれるからです。

7つのカテゴリーに7つ目のダミー変数を加えようとする初心者は大勢います。そうすると回帰ができなくなるので、「ダミー変数の罠」と呼ばれています。切片が含まれている場合、推計するパラメータが1個余分なのです。一般ルールは、カテゴリー数よりも1個少ないダミー変数を使うことです。ですから、四半期データでは3個のダミー変数を、月次データでは11個のダミー変数を、日次データでは6個のダミー変数を、という具合に使います。

各ダミー変数の係数の解釈は、省略したカテゴリーと比較した当該カテゴリーの影響を測るもの、となります。上の例では、月曜日の係数\(d_{1,t}\)は、予測対象変数への日曜日の影響と比べた月曜日の影響を測ることになります。オーストラリアのビール生産量の四半期季節性を取り込むダミー変数の係数の推計値を解釈する例は、以下のようになります。

TSLM()関数は、特殊関数season()を指定すれば、この状況を自動的に取り扱います。

事例: オーストラリアの四半期ビール生産量

7.14に再掲示した、オーストラリア四半期ビール生産量データを思い起こしましょう。

recent_production <- aus_production %>%
  filter(year(Quarter) >= 1992)
recent_production %>%
  autoplot(Beer) +
  labs(y = "メガリットル",
       title = "オーストラリア四半期ビール生産量")
オーストラリア四半期ビール生産量

図 7.14: オーストラリア四半期ビール生産量

将来のビール生産量を予測したいとします。線形トレンドと四半期ダミー変数を持つ回帰モデルを使って、このデータを以下のようにモデル化できます。 \[ y_{t} = \beta_{0} + \beta_{1} t + \beta_{2}d_{2,t} + \beta_3 d_{3,t} + \beta_4 d_{4,t} + \varepsilon_{t} \] ただし、\(t\)が四半期\(i\)なら\(d_{i,t} = 1\)、それ以外なら0です。第1四半期が省略されているので、他の四半期の係数は第1四半期との差を測るものです。

fit_beer <- recent_production %>%
  model(TSLM(Beer ~ trend() + season()))
report(fit_beer)
#> Series: Beer 
#> Model: TSLM 
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -42.90  -7.60  -0.46   7.99  21.79 
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   441.8004     3.7335  118.33  < 2e-16 ***
#> trend()        -0.3403     0.0666   -5.11  2.7e-06 ***
#> season()year2 -34.6597     3.9683   -8.73  9.1e-13 ***
#> season()year3 -17.8216     4.0225   -4.43  3.4e-05 ***
#> season()year4  72.7964     4.0230   18.09  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 12.2 on 69 degrees of freedom
#> Multiple R-squared: 0.924,   Adjusted R-squared: 0.92
#> F-statistic:  211 on 4 and 69 DF, p-value: <2e-16

trend()season()は標準的な関数ではないことに注意してください。それらは、TSLM()のモデル式の中でだけ機能する「特殊」関数です。

平均-0.34メガリットル/四半期の下降トレンドがあります。平均すると、第1四半期生産量と比べて、第2四半期は34.7メガリットル少なく、第3四半期は17.8メガリットル少なく、第4四半期は72.8メガリットル多くなっています。

augment(fit_beer) %>%
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Beer, colour = "実績値")) +
  geom_line(aes(y = .fitted, colour = "当てはめ値")) +
  scale_colour_manual(
    values = c("実績値" = "black", "当てはめ値" = "#D55E00")
  ) +
  labs(y = "メガリットル",
       title = "オーストラリア四半期ビール生産量") +
  guides(colour = guide_legend(title = "系列"))
ビール生産量の実績と予測の時間プロット

図 7.15: ビール生産量の実績と予測の時間プロット

augment(fit_beer) %>%
  ggplot(aes(x = Beer, y = .fitted,
             colour = factor(quarter(Quarter)))) +
  geom_point() +
  labs(y = "当てはめ値", x = "実績値",
       title = "オーストラリア四半期ビール生産量") +
  geom_abline(intercept = 0, slope = 1) +
  guides(colour = guide_legend(title = "四半期"))
ビール生産量の実績と予測の散布図

図 7.16: ビール生産量の実績と予測の散布図

介入変数

予測対象変数に影響があり得る介入をモデル化することは、しばしば必要です。例えば、競合他社の活動、広告費、業界の行動などは全て影響があり得ます。

影響が1期間だけしか続かない場合、「スパイク」変数を使います。これは、介入がある期間では1を、それ以外では0の値を取るダミー変数です。スパイク変数は、外れ値を取り扱うダミー変数と同じです。

中期的や永続的な影響がある介入もあります。介入が水準をシフトさせる(つまり、系列の値が突然かつ永続的に介入時から変化する)なら、「ステップ」変数を使います。ステップ変数は介入以前では0を、介入以降では1の値を取ります。

永続的影響のもう一つの形は、傾きの変化です。その場合、介入時に屈折するトレンド、つまり非線形なトレンドになります。そうした介入は、区分線形トレンドを使って取り扱います。7.7節で議論します。

営業日

営業日の数は月によってかなり違っていることがあり、売上高データに大きな影響を及ぼし得ます。モデル化するには、各月の営業日の数を予測変数として含めることができます。

曜日の数が異なる影響をモデル化できる別のやり方は、以下の予測変数を含めることです。 \[\begin{align*} x_{1} &= \text{月内の月曜の数} \\ x_{2} &= \text{月内の火曜の数} \\ & \vdots \\ x_{7} &= \text{月内の日曜の数} \end{align*}\]

分布ラグ

広告費を予測変数に含めることはしばしば役に立ちます。しかし、広告の影響はキャンペーン終了後も続き得るので、広告費のラグ値も含める必要があります。ですから、以下のような予測変数を使います。 \[\begin{align*} x_{1} &= \text{前月の広告費} \\ x_{2} &= \text{2カ月前の広告費} \\ & \vdots \\ x_{m} &= \text{$m$カ月前の広告費} \end{align*}\]

ラグ数が増えるほど係数は小さくするのが普通です。ただ、これは本書の範囲を超えます。

イースター

イースターがほとんどの祝日と異なるのは、毎年同じ日でない上に、その影響が数日続く点です。このケースでは、祝日になった特定の期間では1、それ以外では0の値を取るダミー変数を使います。

月次データでは、イースターが3月の場合3月に1、イースターが4月の場合4月に1、の値を取るダミー変数を使います。イースターが3月に始まり4月に終わる場合は、ダミー変数をこの2つの月の間で比例的に分割します。

フーリエ系列

季節ダミー変数の替わりに、特に季節周期が長い場合に、使えるのがフーリエ項です。 Jean-Baptiste Fourierは、1700年代生まれのフランスの数学者で、サイン項とコサイン項の系列で、周期を正しく設定すれば、どんな周期の関数も近似できることを示しました。それは、季節パターンにも使えます。

\(m\)を季節周期として、最初のいくつかのフーリエ項は、 \[ x_{1,t} = \sin\left(\textstyle\frac{2\pi t}{m}\right), x_{2,t} = \cos\left(\textstyle\frac{2\pi t}{m}\right), \] \[ x_{3,t} = \sin\left(\textstyle\frac{4\pi t}{m}\right), x_{4,t} = \cos\left(\textstyle\frac{4\pi t}{m}\right), \]

\[ x_{5,t} = \sin\left(\textstyle\frac{6\pi t}{m}\right), x_{6,t} = \cos\left(\textstyle\frac{6\pi t}{m}\right), \] という格好で続きます。季節性が月次なら、これらの予測変数の最初の11個を使って、11個のダミー変数を使ったのと全く同じ予測が得られます。

フーリエ項を使うと、特に\(m\)が大きい場合、ダミー変数を使う場合より少ない数の予測変数で済むのが普通です。ですから、例えば、\(m\approx 52\)な週次データで役に立ちます。季節周期が短いデータ(例えば、四半期データ)では、季節ダミーの替わりにフーリエ項を使う利点はほとんどありません。

これらのフーリエ項はfourier()関数を使って生成できます。例えば、オーストラリアの四半期ビール生産量データは以下のようにモデル化できます。

fourier_beer <- recent_production %>%
  model(TSLM(Beer ~ trend() + fourier(K = 2)))
report(fourier_beer)
#> Series: Beer 
#> Model: TSLM 
#> 
#> Residuals:
#>    Min     1Q Median     3Q    Max 
#> -42.90  -7.60  -0.46   7.99  21.79 
#> 
#> Coefficients:
#>                    Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)        446.8792     2.8732  155.53  < 2e-16 ***
#> trend()             -0.3403     0.0666   -5.11  2.7e-06 ***
#> fourier(K = 2)C1_4   8.9108     2.0112    4.43  3.4e-05 ***
#> fourier(K = 2)S1_4 -53.7281     2.0112  -26.71  < 2e-16 ***
#> fourier(K = 2)C2_4 -13.9896     1.4226   -9.83  9.3e-15 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 12.2 on 69 degrees of freedom
#> Multiple R-squared: 0.924,   Adjusted R-squared: 0.92
#> F-statistic:  211 on 4 and 69 DF, p-value: <2e-16

fourier()の引数Kには、含めるべきサイン項とコサイン項のペア数を指定します。\(m\)を季節周期として、\(K=m/2\)が許される最大です。ここでは最大を使ったので、結果は季節ダミー変数を使って得られたものと同じになっています。

もし最初の2つのフーリエ項(\(x_{1,t}\)\(x_{2,t}\))だけを使うと、季節パターンは単にサインカーブの波に従うことになります。フーリエ項を含む回帰モデルはハーモニック回帰とも言います。最初の2つのフーリエ項の倍音(ハーモニクス)がそれに続くフーリエ項で表されているからです。