16
Help us understand the problem. What are the problem?

posted at

updated at

Statsmodelsによる時系列分析入門

本記事は
image.png

https://www.statsmodels.org/stable/examples/notebooks/generated/autoregressions.html
の記事をDEEPLを用いて翻訳し、適時修正したものです。

StatsmodelsはPythonというプログラミング言語上で動く統計解析ソフトです。statsmodelsのサンプルを動かすにはPCにPythonがインストールされている必要があります。まだインストールされていない方はJupyter notebookのインストールを参照してください。Jupyter notebookはstatsmodelsを動かすのに大変便利です。

勉強会を行います。
2021/03/09(火) 線形回帰・時系列分析入門 無料オンライン勉強会
Python3ではじめる統計学 18回目 初心者大歓迎
https://kinyuzaimu.connpass.com/event/204862/

また、線形回帰に興味のある方はstatsmodelsによる線形回帰 入門も参考にしてみてください。

自己回帰

AutoRegモデルを用いた自己回帰モデルを紹介します。また、AICのような情報基準を最小化するモデルを選択する際に役立つar_select_orderもカバーしています。自己回帰モデルは、次式で与えられます。

image.png

AutoRegの引数は:

  • 決定論的項(トレンド

n: 決定論的な項無し

c: 定数(デフォルト)

ct. 定数と時間トレンド

t: 時間トレンドのみ

  • 季節性のダミー変数(季節性)

Trueはs-1ダミーを含み、sは時系列の期間(例えば、月次の場合は12)である。

  • オーダーメイドの決定論的な項(決定論的項)

  • DeterministicProcessを受け入れる

  • 外生変数 (exog)

  • DataFrameまたは配列の外生変数をモデルに含める

  • 選択されたラグ(lags)の省略

  • lagsが整数の反復可能なものであれば、これらのものだけがモデルに含まれる。

完全な仕様は

image.png

ここで

  • $d_i$ は,mod(t,period)=i の場合に 1 となる季節ダミーである.モデルが定数(c)を含む場合,期間 0 は除外される.

  • t は,最初の観測で 1 から始まる時間トレンド(1,2,...)である。

  • $x_{t,j}$ は外生的回帰変数である.モデルを定義する際には、これらは左側の変数に時間的に配置されていることに注意。

  • $ϵ_t$ はホワイトノイズ過程であると仮定。

つぎの最初のセルは、標準パッケージをインポートし、インラインで表示されるように設定します。

%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import pandas_datareader as pdr
import seaborn as sns
from statsmodels.tsa.ar_model import AutoReg, ar_select_order
from statsmodels.tsa.api import acf, pacf, graphics

つぎのセルは、プロットスタイルを設定し、matplotlib用のpandasの日付変換器を登録し、デフォルトの図形サイズを設定します。

sns.set_style('darkgrid')
pd.plotting.register_matplotlib_converters()
# Default figure size
sns.mpl.rc('figure',figsize=(16, 6))

最初の例では、季節調整されていない米国の住宅着工件数の前月比成長率を使用しています。季節性は、ピークとトラフの規則的なパターンから明らかです。AutoRegを使用する際の警告を避けるために、時系列の頻度を "MS" (month-start)に設定しました。

data = pdr.get_data_fred('HOUSTNSA', '1959-01-01', '2019-06-01')
housing = data.HOUSTNSA.pct_change().dropna()
# Scale by 100 to get percentages
housing = 100 * housing.asfreq('MS')
fig, ax = plt.subplots()
ax = housing.plot(ax=ax)

image.png

まずはAR(3)から始めましょう。これはこのデータを表現する良いモデルではありませんが、APIの基本的な使い方を示しています。

mod = AutoReg(housing, 3, old_names=False)
res = mod.fit()
print(res.summary())

image.png

AutoRegは、OLSと同じ共分散推定器をサポートしています。以下では、Whiteの共分散推定量であるcov_type="HC0 "を使用しています。パラメータ推定値は同じですが、標準誤差に依存する量はすべて変化します。

res = mod.fit(cov_type="HC0")
print(res.summary())

image.png

sel = ar_select_order(housing, 13, old_names=False)
sel.ar_lags
res = sel.model.fit()
print(res.summary())

image.png
image.png

plot_predictは予測を可視化します。ここでは、モデルによって捕捉された文字列の季節性を示す多数の予測を生成しています。

fig = res.plot_predict(720, 840)

image.png

plot_diagnositcsは、モデルがデータの主要な特徴を捉えていることを示します。

fig = plt.figure(figsize=(16,9))
fig = res.plot_diagnostics(fig=fig, lags=30)

image.png

季節性ダミー変数

AutoRegは、季節性をモデル化する別の方法である季節ダミーをサポートしています。ダミーを含めることで、AR(2)だけのダイナミクスに短縮されます。

AutoRegは、季節性をモデル化する季節ダミーをサポートしています。ダミーを含めることで、AR(2)だけのダイナミクスに短縮されます。

sel = ar_select_order(housing, 13, seasonal=True, old_names=False)
sel.ar_lags
res = sel.model.fit()
print(res.summary())
#
AutoReg Model Results                             
==============================================================================
Dep. Variable:               HOUSTNSA   No. Observations:                  725
Model:               Seas. AutoReg(2)   Log Likelihood               -2652.556
Method:               Conditional MLE   S.D. of innovations              9.487
Date:                Tue, 23 Feb 2021   AIC                              4.541
Time:                        14:04:37   BIC                              4.636
Sample:                    04-01-1959   HQIC                             4.578
                         - 06-01-2019                                         
===============================================================================
                  coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------------------------------------------
const           1.2726      1.373      0.927      0.354      -1.418       3.963
s(2,12)        32.6477      1.824     17.901      0.000      29.073      36.222
s(3,12)        23.0685      2.435      9.472      0.000      18.295      27.842
s(4,12)        10.7267      2.693      3.983      0.000       5.449      16.005
s(5,12)         1.6792      2.100      0.799      0.424      -2.437       5.796
s(6,12)        -4.4229      1.896     -2.333      0.020      -8.138      -0.707
s(7,12)        -4.2113      1.824     -2.309      0.021      -7.786      -0.636
s(8,12)        -6.4124      1.791     -3.581      0.000      -9.922      -2.902
s(9,12)         0.1095      1.800      0.061      0.952      -3.419       3.638
s(10,12)      -16.7511      1.814     -9.234      0.000     -20.307     -13.196
s(11,12)      -20.7023      1.862    -11.117      0.000     -24.352     -17.053
s(12,12)      -11.9554      1.778     -6.724      0.000     -15.440      -8.470
HOUSTNSA.L1    -0.2953      0.037     -7.994      0.000      -0.368      -0.223
HOUSTNSA.L2    -0.1148      0.037     -3.107      0.002      -0.187      -0.042
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1           -1.2862           -2.6564j            2.9514           -0.3218
AR.2           -1.2862           +2.6564j            2.9514            0.3218
-----------------------------------------------------------------------------

季節的なダミーは、10年後の未来のすべての期間に自明な季節的要素を持たない予測を見れば一目瞭然である。

fig = res.plot_predict(720, 840)

image.png

fig = plt.figure(figsize=(16,9))
fig = res.plot_diagnostics(lags=30, fig=fig)

image.png

Seasonal Dynamics

AutoRegはOLSを使用してパラメータを推定するため、季節成分を直接サポートしていませんが、季節ARの制約を課さないオーバーパラメトリックな季節ARを使用して季節のダイナミクスを捉えることができます。

yoy_housing = data.HOUSTNSA.pct_change(12).resample("MS").last().dropna()
_, ax = plt.subplots()
ax = yoy_housing.plot(ax=ax)

image.png
最大ラグのみを選択する単純な手法を用いてモデルを選択することから始めます。すべての下位ラグは自動的に含まれます。これは、モデルが短期AR(1)成分と季節的AR(1)成分の両方を持つ季節的ARの次を可能にするので、チェックすべき最大ラグは13に設定されています。
image.png

を展開したときに AutoRegは構造を強制しませんが、ネスティングモデルを推定することができます。
image.png
13 のラグがすべて選択されていることがわかります。

sel = ar_select_order(yoy_housing, 13, old_names=False)
sel.ar_lags
#
array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13])

13本のラグがすべて必要とされることはなさそうです。glob=Trueを設定して、最大13のラグを含む213のモデルをすべて検索することができます。

ここでは、最初の3つが選択され、7番目が選択され、最後に12番目と13番目が選択されていることがわかります。これは、表面的には上述の構造に似ています。

モデルをフィッティングした後、この仕様がデータのダイナミクスを捉えるのに十分であるように見えることを示す診断プロットを見てみましょう。

sel = ar_select_order(yoy_housing, 13, glob=True, old_names=False)
sel.ar_lags
res = sel.model.fit()
print(res.summary())

#
AutoReg Model Results                             
==============================================================================
Dep. Variable:               HOUSTNSA   No. Observations:                  714
Model:             Restr. AutoReg(13)   Log Likelihood                 589.177
Method:               Conditional MLE   S.D. of innovations              0.104
Date:                Tue, 23 Feb 2021   AIC                             -4.496
Time:                        14:06:37   BIC                             -4.444
Sample:                    02-01-1961   HQIC                            -4.476
                         - 06-01-2019                                         
================================================================================
                   coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------
const            0.0035      0.004      0.875      0.382      -0.004       0.011
HOUSTNSA.L1      0.5640      0.035     16.167      0.000       0.496       0.632
HOUSTNSA.L2      0.2347      0.038      6.238      0.000       0.161       0.308
HOUSTNSA.L3      0.2051      0.037      5.560      0.000       0.133       0.277
HOUSTNSA.L7     -0.0903      0.030     -2.976      0.003      -0.150      -0.031
HOUSTNSA.L12    -0.3791      0.034    -11.075      0.000      -0.446      -0.312
HOUSTNSA.L13     0.3354      0.033     10.254      0.000       0.271       0.400
                                    Roots                                     
==============================================================================
                   Real          Imaginary           Modulus         Frequency
------------------------------------------------------------------------------
AR.1            -1.0309           -0.2682j            1.0652           -0.4595
AR.2            -1.0309           +0.2682j            1.0652            0.4595
AR.3            -0.7454           -0.7417j            1.0515           -0.3754
AR.4            -0.7454           +0.7417j            1.0515            0.3754
AR.5            -0.3172           -1.0221j            1.0702           -0.2979
AR.6            -0.3172           +1.0221j            1.0702            0.2979
AR.7             0.2419           -1.0573j            1.0846           -0.2142
AR.8             0.2419           +1.0573j            1.0846            0.2142
AR.9             0.7840           -0.8303j            1.1420           -0.1296
AR.10            0.7840           +0.8303j            1.1420            0.1296
AR.11            1.0730           -0.2386j            1.0992           -0.0348
AR.12            1.0730           +0.2386j            1.0992            0.0348
AR.13            1.1193           -0.0000j            1.1193           -0.0000
------------------------------------------------------------------------------
fig = plt.figure(figsize=(16,9))
fig = res.plot_diagnostics(fig=fig, lags=30)

image.png

季節的なダミーも含めることができます。モデルは前年比変化を使用しているので、これらはすべて有意ではありません。

sel = ar_select_order(yoy_housing, 13, glob=True, seasonal=True, old_names=False)
sel.ar_lags
res = sel.model.fit()
print(res.summary())

#
AutoReg Model Results                                
====================================================================================
Dep. Variable:                     HOUSTNSA   No. Observations:                  714
Model:             Restr. Seas. AutoReg(13)   Log Likelihood                 590.875
Method:                     Conditional MLE   S.D. of innovations              0.104
Date:                      Tue, 23 Feb 2021   AIC                             -4.469
Time:                              14:07:18   BIC                             -4.346
Sample:                          02-01-1961   HQIC                            -4.422
                               - 06-01-2019                                         
================================================================================
                   coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------
const            0.0167      0.014      1.215      0.224      -0.010       0.044
s(2,12)         -0.0179      0.019     -0.931      0.352      -0.056       0.020
s(3,12)         -0.0121      0.019     -0.630      0.528      -0.050       0.026
s(4,12)         -0.0210      0.019     -1.089      0.276      -0.059       0.017
s(5,12)         -0.0223      0.019     -1.157      0.247      -0.060       0.015
s(6,12)         -0.0224      0.019     -1.160      0.246      -0.060       0.015
s(7,12)         -0.0212      0.019     -1.096      0.273      -0.059       0.017
s(8,12)         -0.0101      0.019     -0.520      0.603      -0.048       0.028
s(9,12)         -0.0095      0.019     -0.491      0.623      -0.047       0.028
s(10,12)        -0.0049      0.019     -0.252      0.801      -0.043       0.033
s(11,12)        -0.0084      0.019     -0.435      0.664      -0.046       0.030
s(12,12)        -0.0077      0.019     -0.400      0.689      -0.046       0.030
HOUSTNSA.L1      0.5630      0.035     16.160      0.000       0.495       0.631
HOUSTNSA.L2      0.2347      0.038      6.248      0.000       0.161       0.308
HOUSTNSA.L3      0.2075      0.037      5.634      0.000       0.135       0.280
HOUSTNSA.L7     -0.0916      0.030     -3.013      0.003      -0.151      -0.032
HOUSTNSA.L12    -0.3810      0.034    -11.149      0.000      -0.448      -0.314
HOUSTNSA.L13     0.3373      0.033     10.327      0.000       0.273       0.401
                                    Roots                                     
==============================================================================
                   Real          Imaginary           Modulus         Frequency
------------------------------------------------------------------------------
AR.1            -1.0305           -0.2681j            1.0648           -0.4595
AR.2            -1.0305           +0.2681j            1.0648            0.4595
AR.3            -0.7447           -0.7414j            1.0509           -0.3754
AR.4            -0.7447           +0.7414j            1.0509            0.3754
AR.5            -0.3172           -1.0215j            1.0696           -0.2979
AR.6            -0.3172           +1.0215j            1.0696            0.2979
AR.7             0.2416           -1.0568j            1.0841           -0.2142
AR.8             0.2416           +1.0568j            1.0841            0.2142
AR.9             0.7837           -0.8304j            1.1418           -0.1296
AR.10            0.7837           +0.8304j            1.1418            0.1296
AR.11            1.0724           -0.2383j            1.0986           -0.0348
AR.12            1.0724           +0.2383j            1.0986            0.0348
AR.13            1.1192           -0.0000j            1.1192           -0.0000
------------------------------------------------------------------------------

工業生産

工業生産指数データを用いて予測を検討します。

data = pdr.get_data_fred('INDPRO', '1959-01-01', '2019-06-01')
ind_prod = data.INDPRO.pct_change(12).dropna().asfreq('MS')
_, ax = plt.subplots(figsize=(16,9))
ind_prod.plot(ax=ax)

image.png
最大12のラグを使用するモデルを選択することから始めます。AR(13) は、多くの係数が取るに足らないものであっても、BIC 基準を最小化します。

sel = ar_select_order(ind_prod, 13, 'bic', old_names=False)
res = sel.model.fit()
print(res.summary())
"
AutoReg Model Results                             
==============================================================================
Dep. Variable:                 INDPRO   No. Observations:                  714
Model:                    AutoReg(13)   Log Likelihood                2318.848
Method:               Conditional MLE   S.D. of innovations              0.009
Date:                Tue, 23 Feb 2021   AIC                             -9.411
Time:                        14:08:08   BIC                             -9.313
Sample:                    02-01-1961   HQIC                            -9.373
                         - 06-01-2019                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.0012      0.000      2.896      0.004       0.000       0.002
INDPRO.L1      1.1504      0.035     32.940      0.000       1.082       1.219
INDPRO.L2     -0.0747      0.053     -1.407      0.159      -0.179       0.029
INDPRO.L3      0.0043      0.053      0.081      0.935      -0.099       0.107
INDPRO.L4      0.0027      0.053      0.052      0.959      -0.100       0.106
INDPRO.L5     -0.1383      0.052     -2.642      0.008      -0.241      -0.036
INDPRO.L6      0.0085      0.052      0.163      0.871      -0.094       0.111
INDPRO.L7      0.0375      0.052      0.720      0.471      -0.065       0.139
INDPRO.L8     -0.0235      0.052     -0.453      0.651      -0.125       0.078
INDPRO.L9      0.0945      0.052      1.824      0.068      -0.007       0.196
INDPRO.L10    -0.0844      0.052     -1.627      0.104      -0.186       0.017
INDPRO.L11     0.0024      0.052      0.047      0.962      -0.099       0.104
INDPRO.L12    -0.3809      0.052     -7.367      0.000      -0.482      -0.280
INDPRO.L13     0.3589      0.033     10.916      0.000       0.294       0.423
                                    Roots                                     
==============================================================================
                   Real          Imaginary           Modulus         Frequency
------------------------------------------------------------------------------
AR.1            -1.0411           -0.2889j            1.0804           -0.4569
AR.2            -1.0411           +0.2889j            1.0804            0.4569
AR.3            -0.7774           -0.8064j            1.1201           -0.3721
AR.4            -0.7774           +0.8064j            1.1201            0.3721
AR.5            -0.2760           -1.0530j            1.0886           -0.2908
AR.6            -0.2760           +1.0530j            1.0886            0.2908
AR.7             0.2719           -1.0515j            1.0861           -0.2097
AR.8             0.2719           +1.0515j            1.0861            0.2097
AR.9             0.8019           -0.7297j            1.0842           -0.1175
AR.10            0.8019           +0.7297j            1.0842            0.1175
AR.11            1.0224           -0.2214j            1.0461           -0.0339
AR.12            1.0224           +0.2214j            1.0461            0.0339
AR.13            1.0578           -0.0000j            1.0578           -0.0000
------------------------------------------------------------------------------

また、短いラグを必要とせずに、必要に応じてより長いラグを入力できるグローバル検索を使用することもできます。ここでは、多くのラグが削除されているのがわかります。このモデルは、データに季節性があるかもしれないことを示しています。

sel = ar_select_order(ind_prod, 13, 'bic', glob=True, old_names=False)
sel.ar_lags
res_glob = sel.model.fit()
print(res.summary())
#
 AutoReg Model Results                             
==============================================================================
Dep. Variable:                 INDPRO   No. Observations:                  714
Model:                    AutoReg(13)   Log Likelihood                2318.848
Method:               Conditional MLE   S.D. of innovations              0.009
Date:                Tue, 23 Feb 2021   AIC                             -9.411
Time:                        14:08:26   BIC                             -9.313
Sample:                    02-01-1961   HQIC                            -9.373
                         - 06-01-2019                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.0012      0.000      2.896      0.004       0.000       0.002
INDPRO.L1      1.1504      0.035     32.940      0.000       1.082       1.219
INDPRO.L2     -0.0747      0.053     -1.407      0.159      -0.179       0.029
INDPRO.L3      0.0043      0.053      0.081      0.935      -0.099       0.107
INDPRO.L4      0.0027      0.053      0.052      0.959      -0.100       0.106
INDPRO.L5     -0.1383      0.052     -2.642      0.008      -0.241      -0.036
INDPRO.L6      0.0085      0.052      0.163      0.871      -0.094       0.111
INDPRO.L7      0.0375      0.052      0.720      0.471      -0.065       0.139
INDPRO.L8     -0.0235      0.052     -0.453      0.651      -0.125       0.078
INDPRO.L9      0.0945      0.052      1.824      0.068      -0.007       0.196
INDPRO.L10    -0.0844      0.052     -1.627      0.104      -0.186       0.017
INDPRO.L11     0.0024      0.052      0.047      0.962      -0.099       0.104
INDPRO.L12    -0.3809      0.052     -7.367      0.000      -0.482      -0.280
INDPRO.L13     0.3589      0.033     10.916      0.000       0.294       0.423
                                    Roots                                     
==============================================================================
                   Real          Imaginary           Modulus         Frequency
------------------------------------------------------------------------------
AR.1            -1.0411           -0.2889j            1.0804           -0.4569
AR.2            -1.0411           +0.2889j            1.0804            0.4569
AR.3            -0.7774           -0.8064j            1.1201           -0.3721
AR.4            -0.7774           +0.8064j            1.1201            0.3721
AR.5            -0.2760           -1.0530j            1.0886           -0.2908
AR.6            -0.2760           +1.0530j            1.0886            0.2908
AR.7             0.2719           -1.0515j            1.0861           -0.2097
AR.8             0.2719           +1.0515j            1.0861            0.2097
AR.9             0.8019           -0.7297j            1.0842           -0.1175
AR.10            0.8019           +0.7297j            1.0842            0.1175
AR.11            1.0224           -0.2214j            1.0461           -0.0339
AR.12            1.0224           +0.2214j            1.0461            0.0339
AR.13            1.0578           -0.0000j            1.0578           -0.0000
------------------------------------------------------------------------------

plot_predictは、信頼区間とともに予測プロットを作成します。ここでは、最後の観測から始まり、18ヶ月間の予測を作成しています。

ind_prod.shape
#
(714,)
fig = res_glob.plot_predict(start=714, end=732)

image.png
フルモデルと制限モデルからの予測は非常によく似ています。また、非常に異なるダイナミクスを持つAR(5)も含みます。

res_ar5 = AutoReg(ind_prod, 5, old_names=False).fit()
predictions = pd.DataFrame({"AR(5)": res_ar5.predict(start=714, end=726),
                            "AR(13)": res.predict(start=714, end=726),
                            "Restr. AR(13)": res_glob.predict(start=714, end=726)})
_, ax = plt.subplots()
ax = predictions.plot(ax=ax)

image.png
診断は、モデルがデータのダイナミクスの大部分を捉えていることを示しています。ACFは季節周波数でのパターンを示しているので、より完全な季節モデル(SARIMAX)が必要かもしれません。

fig = plt.figure(figsize=(16,9))
fig = res_glob.plot_diagnostics(fig=fig, lags=30)

image.png

予測

予測は、結果のインスタンスから predict メソッドを使用して生成されます。デフォルトでは、静的な予測が生成され、これはワンステップ予測です。マルチステップ予測を作成するには、dynamic=Trueを使用する必要があります。

この次のセルでは、サンプルの最後の24期間の12段階の予測を作成します。これにはループが必要です。

注意:これらは、我々が予測しているデータがパラメータを推定するために使用されたので、技術的にはインサンプルです。OOS予測を作成するには、2つのモデルが必要です。1つ目は、OOS期間を除外しなければなりません。2つ目は、OOS期間を除外した短いサンプル・モデルのパラメータを用いて、フル・サンプル・モデルからの予測法を使用します。

import numpy as np
start = ind_prod.index[-24]
forecast_index = pd.date_range(start, freq=ind_prod.index.freq, periods=36)
cols = ['-'.join(str(val) for val in (idx.year, idx.month)) for idx in forecast_index]
forecasts = pd.DataFrame(index=forecast_index,columns=cols)
for i in range(1, 24):
    fcast = res_glob.predict(start=forecast_index[i], end=forecast_index[i+12], dynamic=True)
    forecasts.loc[fcast.index, cols[i]] = fcast
_, ax = plt.subplots(figsize=(16, 10))
ind_prod.iloc[-24:].plot(ax=ax, color="black", linestyle="--")
ax = forecasts.plot(ax=ax)

image.png

Comparing to SARIMAX

from statsmodels.tsa.api import SARIMAX

sarimax_mod = SARIMAX(ind_prod, order=((1,5,12,13),0, 0), trend='c')
sarimax_res = sarimax_mod.fit()
print(sarimax_res.summary())

sarimax_params = sarimax_res.params.iloc[:-1].copy()
sarimax_params.index = res_glob.params.index
params = pd.concat([res_glob.params, sarimax_params], axis=1, sort=False)
params.columns = ["AutoReg", "SARIMAX"]
params

SARIMAXとの比較

SARIMAXは、eXogenous regressorsモデルを用いた季節自己回帰統合移動平均の実装です。SARIMAXは以下をサポートしています。

  • 季節的および非季節的なARおよびMA成分の指定

  • 外生変数を含める

  • カルマンフィルタを用いた完全最尤推定

このモデルはAutoRegよりも特徴が豊富です。SARIMAXとは異なり、AutoRegはOLSを用いてパラメータを推定します。これはより高速であり、問題は大域的に凸であるため、局所的なミニマの問題はありません。AR(P)モデルを比較する際に、閉形式推定器とその性能がSARIMAXに対するAutoRegの主な利点となります。AutoRegは季節ダミーもサポートしており、ユーザーが外生的回帰子としてそれらを含める場合、SARIMAXと一緒に使用することができます。

from statsmodels.tsa.api import SARIMAX

sarimax_mod = SARIMAX(ind_prod, order=((1,5,12,13),0, 0), trend='c')
sarimax_res = sarimax_mod.fit()
print(sarimax_res.summary())
#
SARIMAX Results                                     
=========================================================================================
Dep. Variable:                            INDPRO   No. Observations:                  714
Model:             SARIMAX([1, 5, 12, 13], 0, 0)   Log Likelihood                2302.146
Date:                           Tue, 23 Feb 2021   AIC                          -4592.293
Time:                                   14:10:42   BIC                          -4564.867
Sample:                               01-01-1960   HQIC                         -4581.701
                                    - 06-01-2019                                         
Covariance Type:                             opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      0.0011      0.000      2.627      0.009       0.000       0.002
ar.L1          1.0798      0.010    106.388      0.000       1.060       1.100
ar.L5         -0.0849      0.011     -7.467      0.000      -0.107      -0.063
ar.L12        -0.4392      0.026    -16.935      0.000      -0.490      -0.388
ar.L13         0.4032      0.025     15.831      0.000       0.353       0.453
sigma2      9.179e-05   3.13e-06     29.366      0.000    8.57e-05    9.79e-05
===================================================================================
Ljung-Box (L1) (Q):                  19.65   Jarque-Bera (JB):               928.85
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               0.38   Skew:                            -0.63
Prob(H) (two-sided):                  0.00   Kurtosis:                         8.44
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
sarimax_params = sarimax_res.params.iloc[:-1].copy()
sarimax_params.index = res_glob.params.index
params = pd.concat([res_glob.params, sarimax_params], axis=1, sort=False)
params.columns = ["AutoReg", "SARIMAX"]
params
#
AutoReg SARIMAX
const   0.001292    0.001140
INDPRO.L1   1.088215    1.079829
INDPRO.L5   -0.105668   -0.084904
INDPRO.L12  -0.385171   -0.439175
INDPRO.L13  0.358738    0.403241

カスタム決定論的プロセス

決定論的パラメータは、自由に決定論的プロセスを作成ことを可能にします。これにより、より複雑な決定論的項を構築することができます。例えば、2つの期間を持つ季節成分を含むものや、次の例が示すように、季節ダミーではなくフーリエ級数を使用するものなどがあります。

from statsmodels.tsa.deterministic import DeterministicProcess
dp = DeterministicProcess(housing.index, constant=True, period=12, fourier=2)
mod = AutoReg(housing,2, trend="n",seasonal=False, deterministic=dp)
res = mod.fit()
print(res.summary())
#
 AutoReg Model Results                             
==============================================================================
Dep. Variable:               HOUSTNSA   No. Observations:                  725
Model:                     AutoReg(2)   Log Likelihood               -2716.505
Method:               Conditional MLE   S.D. of innovations             10.364
Date:                Tue, 23 Feb 2021   AIC                              4.699
Time:                        14:11:45   BIC                              4.750
Sample:                    04-01-1959   HQIC                             4.718
                         - 06-01-2019                                         
===============================================================================
                  coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------------------------------------------
const           1.7550      0.391      4.485      0.000       0.988       2.522
sin(1,12)      16.7443      0.860     19.478      0.000      15.059      18.429
cos(1,12)       4.9409      0.588      8.409      0.000       3.789       6.093
sin(2,12)      12.9364      0.619     20.889      0.000      11.723      14.150
cos(2,12)      -0.4738      0.754     -0.628      0.530      -1.952       1.004
HOUSTNSA.L1    -0.3905      0.037    -10.664      0.000      -0.462      -0.319
HOUSTNSA.L2    -0.1746      0.037     -4.769      0.000      -0.246      -0.103
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1           -1.1182           -2.1159j            2.3932           -0.3274
AR.2           -1.1182           +2.1159j            2.3932            0.3274
-----------------------------------------------------------------------------
fig = res.plot_predict(720, 840)
fig = res.plot_predict(720, 840)

image.png

確定的トレンド

基本的な使い方

基本的な構成は DeterministicProcess を通して直接構築することができます。これらは、定数、任意の順序の時間トレンド、季節成分またはフーリエ成分で構成されます。

プロセスは、フルサンプル(またはインサンプル)のインデックスを必要とします。

まず、定数、線形時間トレンド、5周期の季節項を持つ決定論的プロセスを初期化します。in_sampleメソッドは、インデックスに一致する値の完全集合を返します。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

plt.rc("figure", figsize=(16, 9))
plt.rc("font", size=16)
from statsmodels.tsa.deterministic import DeterministicProcess

index = pd.RangeIndex(0, 100)
det_proc = DeterministicProcess(
    index, constant=True, order=1, seasonal=True, period=5
)
det_proc.in_sample()
#
const   trend   s(2,5)  s(3,5)  s(4,5)  s(5,5)
0   1.0 1.0 0.0 0.0 0.0 0.0
1   1.0 2.0 1.0 0.0 0.0 0.0
2   1.0 3.0 0.0 1.0 0.0 0.0
3   1.0 4.0 0.0 0.0 1.0 0.0
4   1.0 5.0 0.0 0.0 0.0 1.0
... ... ... ... ... ... ...
95  1.0 96.0    0.0 0.0 0.0 0.0
96  1.0 97.0    1.0 0.0 0.0 0.0
97  1.0 98.0    0.0 1.0 0.0 0.0
98  1.0 99.0    0.0 0.0 1.0 0.0
99  1.0 100.0   0.0 0.0 0.0 1.0
100 rows × 6 columns

out_of_sample は、インサンプル終了後の次のステップの値を返します。

det_proc.out_of_sample(15)
#
const   trend   s(2,5)  s(3,5)  s(4,5)  s(5,5)
100 1.0 101.0   0.0 0.0 0.0 0.0
101 1.0 102.0   1.0 0.0 0.0 0.0
102 1.0 103.0   0.0 1.0 0.0 0.0
103 1.0 104.0   0.0 0.0 1.0 0.0
104 1.0 105.0   0.0 0.0 0.0 1.0
105 1.0 106.0   0.0 0.0 0.0 0.0
106 1.0 107.0   1.0 0.0 0.0 0.0
107 1.0 108.0   0.0 1.0 0.0 0.0
108 1.0 109.0   0.0 0.0 1.0 0.0
109 1.0 110.0   0.0 0.0 0.0 1.0
110 1.0 111.0   0.0 0.0 0.0 0.0
111 1.0 112.0   1.0 0.0 0.0 0.0
112 1.0 113.0   0.0 1.0 0.0 0.0
113 1.0 114.0   0.0 0.0 1.0 0.0
114 1.0 115.0   0.0 0.0 0.0 1.0

range(start, stop) は、サンプル内およびサンプル外を含む任意の範囲の決定論的項を生成するために使用できます。

注意事項

インデックスがpandas DatetimeIndexまたはPeriodIndexである場合、startとstopは日付的なもの(文字列、例えば "2020-06-01 "やTimestamp)でも整数でも構いません。

stopは常に範囲に含まれます。これはあまりPythonicなものではありませんが、statsmodelsとPandasの両方がデータライクなスライスを扱うときにstopを含むので必要です。

det_proc.range(190, 210)
#
const   trend   s(2,5)  s(3,5)  s(4,5)  s(5,5)
190 1.0 191.0   0.0 0.0 0.0 0.0
191 1.0 192.0   1.0 0.0 0.0 0.0
192 1.0 193.0   0.0 1.0 0.0 0.0
193 1.0 194.0   0.0 0.0 1.0 0.0
194 1.0 195.0   0.0 0.0 0.0 1.0
195 1.0 196.0   0.0 0.0 0.0 0.0
196 1.0 197.0   1.0 0.0 0.0 0.0
197 1.0 198.0   0.0 1.0 0.0 0.0
198 1.0 199.0   0.0 0.0 1.0 0.0
199 1.0 200.0   0.0 0.0 0.0 1.0
200 1.0 201.0   0.0 0.0 0.0 0.0
201 1.0 202.0   1.0 0.0 0.0 0.0
202 1.0 203.0   0.0 1.0 0.0 0.0
203 1.0 204.0   0.0 0.0 1.0 0.0
204 1.0 205.0   0.0 0.0 0.0 1.0
205 1.0 206.0   0.0 0.0 0.0 0.0
206 1.0 207.0   1.0 0.0 0.0 0.0
207 1.0 208.0   0.0 1.0 0.0 0.0
208 1.0 209.0   0.0 0.0 1.0 0.0
209 1.0 210.0   0.0 0.0 0.0 1.0
210 1.0 211.0   0.0 0.0 0.0 0.0

日付のインデックス

次に、PeriodIndexを使って同じ手順を示します。

index = pd.period_range("2020-03-01", freq="M", periods=60)
det_proc = DeterministicProcess(index, constant=True, fourier=2)
det_proc.in_sample().head(12)
#
const   sin(1,12)   cos(1,12)   sin(2,12)   cos(2,12)
2020-03 1.0 0.000000e+00    1.000000e+00    0.000000e+00    1.0
2020-04 1.0 5.000000e-01    8.660254e-01    8.660254e-01    0.5
2020-05 1.0 8.660254e-01    5.000000e-01    8.660254e-01    -0.5
2020-06 1.0 1.000000e+00    6.123234e-17    1.224647e-16    -1.0
2020-07 1.0 8.660254e-01    -5.000000e-01   -8.660254e-01   -0.5
2020-08 1.0 5.000000e-01    -8.660254e-01   -8.660254e-01   0.5
2020-09 1.0 1.224647e-16    -1.000000e+00   -2.449294e-16   1.0
2020-10 1.0 -5.000000e-01   -8.660254e-01   8.660254e-01    0.5
2020-11 1.0 -8.660254e-01   -5.000000e-01   8.660254e-01    -0.5
2020-12 1.0 -1.000000e+00   -1.836970e-16   3.673940e-16    -1.0
2021-01 1.0 -8.660254e-01   5.000000e-01    -8.660254e-01   -0.5
2021-02 1.0 -5.000000e-01   8.660254e-01    -8.660254e-01   0.5
det_proc.out_of_sample(12)
#
const   sin(1,12)   cos(1,12)   sin(2,12)   cos(2,12)
2025-03 1.0 -1.224647e-15   1.000000e+00    -2.449294e-15   1.0
2025-04 1.0 5.000000e-01    8.660254e-01    8.660254e-01    0.5
2025-05 1.0 8.660254e-01    5.000000e-01    8.660254e-01    -0.5
2025-06 1.0 1.000000e+00    -4.904777e-16   -9.809554e-16   -1.0
2025-07 1.0 8.660254e-01    -5.000000e-01   -8.660254e-01   -0.5
2025-08 1.0 5.000000e-01    -8.660254e-01   -8.660254e-01   0.5
2025-09 1.0 4.899825e-15    -1.000000e+00   -9.799650e-15   1.0
2025-10 1.0 -5.000000e-01   -8.660254e-01   8.660254e-01    0.5
2025-11 1.0 -8.660254e-01   -5.000000e-01   8.660254e-01    -0.5
2025-12 1.0 -1.000000e+00   -3.184701e-15   6.369401e-15    -1.0
2026-01 1.0 -8.660254e-01   5.000000e-01    -8.660254e-01   -0.5
2026-02 1.0 -5.000000e-01   8.660254e-01    -8.660254e-01   0.5

rangeは、通常は文字列として与えられる日付の引数を受け付けます。

det_proc.range("2025-01", "2026-01")
#
const   sin(1,12)   cos(1,12)   sin(2,12)   cos(2,12)
2025-01 1.0 -8.660254e-01   5.000000e-01    -8.660254e-01   -0.5
2025-02 1.0 -5.000000e-01   8.660254e-01    -8.660254e-01   0.5
2025-03 1.0 -1.224647e-15   1.000000e+00    -2.449294e-15   1.0
2025-04 1.0 5.000000e-01    8.660254e-01    8.660254e-01    0.5
2025-05 1.0 8.660254e-01    5.000000e-01    8.660254e-01    -0.5
2025-06 1.0 1.000000e+00    -4.904777e-16   -9.809554e-16   -1.0
2025-07 1.0 8.660254e-01    -5.000000e-01   -8.660254e-01   -0.5
2025-08 1.0 5.000000e-01    -8.660254e-01   -8.660254e-01   0.5
2025-09 1.0 4.899825e-15    -1.000000e+00   -9.799650e-15   1.0
2025-10 1.0 -5.000000e-01   -8.660254e-01   8.660254e-01    0.5
2025-11 1.0 -8.660254e-01   -5.000000e-01   8.660254e-01    -0.5
2025-12 1.0 -1.000000e+00   -3.184701e-15   6.369401e-15    -1.0
2026-01 1.0 -8.660254e-01   5.000000e-01    -8.660254e-01   -0.5

これは、整数値58と整数値70を使用してもできます。

det_proc.range(58, 70)
#
const   sin(1,12)   cos(1,12)   sin(2,12)   cos(2,12)
2025-01 1.0 -8.660254e-01   5.000000e-01    -8.660254e-01   -0.5
2025-02 1.0 -5.000000e-01   8.660254e-01    -8.660254e-01   0.5
2025-03 1.0 -1.224647e-15   1.000000e+00    -2.449294e-15   1.0
2025-04 1.0 5.000000e-01    8.660254e-01    8.660254e-01    0.5
2025-05 1.0 8.660254e-01    5.000000e-01    8.660254e-01    -0.5
2025-06 1.0 1.000000e+00    -4.904777e-16   -9.809554e-16   -1.0
2025-07 1.0 8.660254e-01    -5.000000e-01   -8.660254e-01   -0.5
2025-08 1.0 5.000000e-01    -8.660254e-01   -8.660254e-01   0.5
2025-09 1.0 4.899825e-15    -1.000000e+00   -9.799650e-15   1.0
2025-10 1.0 -5.000000e-01   -8.660254e-01   8.660254e-01    0.5
2025-11 1.0 -8.660254e-01   -5.000000e-01   8.660254e-01    -0.5
2025-12 1.0 -1.000000e+00   -3.184701e-15   6.369401e-15    -1.0
2026-01 1.0 -8.660254e-01   5.000000e-01    -8.660254e-01   -0.5

コントラクタの作成

コンストラクタで直接サポートされていない決定論的プロセスは、 additional_terms を用いて DetermisticTerm のリストを受け取ることで作成できます。ここでは、2つの季節成分を持つ決定論的プロセスを作成しています: 5日間の期間を持つ曜日と、365.25日の期間を持つフーリエ成分で捕捉された年次です。

from statsmodels.tsa.deterministic import Fourier, Seasonality, TimeTrend

index = pd.period_range("2020-03-01", freq="D", periods=2 * 365)
tt = TimeTrend(constant=True)
four = Fourier(period=365.25, order=2)
seas = Seasonality(period=7)
det_proc = DeterministicProcess(index, additional_terms=[tt, seas, four])
det_proc.in_sample().head(28)
#
const   s(2,7)  s(3,7)  s(4,7)  s(5,7)  s(6,7)  s(7,7)  sin(1,365.25)   cos(1,365.25)   sin(2,365.25)   cos(2,365.25)
2020-03-01  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.000000    1.000000    0.000000    1.000000
2020-03-02  1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.017202    0.999852    0.034398    0.999408
2020-03-03  1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.034398    0.999408    0.068755    0.997634
2020-03-04  1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.051584    0.998669    0.103031    0.994678
2020-03-05  1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.068755    0.997634    0.137185    0.990545
2020-03-06  1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.085906    0.996303    0.171177    0.985240
2020-03-07  1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.103031    0.994678    0.204966    0.978769
2020-03-08  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.120126    0.992759    0.238513    0.971139
2020-03-09  1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.137185    0.990545    0.271777    0.962360
2020-03-10  1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.154204    0.988039    0.304719    0.952442
2020-03-11  1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.171177    0.985240    0.337301    0.941397
2020-03-12  1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.188099    0.982150    0.369484    0.929237
2020-03-13  1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.204966    0.978769    0.401229    0.915978
2020-03-14  1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.221772    0.975099    0.432499    0.901634
2020-03-15  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.238513    0.971139    0.463258    0.886224
2020-03-16  1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.255182    0.966893    0.493468    0.869764
2020-03-17  1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.271777    0.962360    0.523094    0.852275
2020-03-18  1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.288291    0.957543    0.552101    0.833777
2020-03-19  1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.304719    0.952442    0.580455    0.814292
2020-03-20  1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.321058    0.947060    0.608121    0.793844
2020-03-21  1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.337301    0.941397    0.635068    0.772456
2020-03-22  1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.353445    0.935455    0.661263    0.750154
2020-03-23  1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.369484    0.929237    0.686676    0.726964
2020-03-24  1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.385413    0.922744    0.711276    0.702913
2020-03-25  1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.401229    0.915978    0.735034    0.678031
2020-03-26  1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.416926    0.908940    0.757922    0.652346
2020-03-27  1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.432499    0.901634    0.779913    0.625889
2020-03-28  1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.447945    0.894061    0.800980    0.598691

決定論的プロセスの作成

DetermisticTerm 抽象基底クラスは、ユーザーが自由に決定論的な項を書くのを助けるためにサブクラス化されています。次に2つの例を示します。1つ目は、一定の期間の後にブレークするブレークタイムトレンドです。2つ目は、決定論的なプロセスではない外生データを、あたかも決定論的であるかのように扱うことを可能にする、"トリック "な決定論的項です。これにより、予測に必要な項の生成を単純化することができます。

これらは、カスタム用語の構築を可能にすることを意図し、入力の妥当性の確認の点で間違いなく良い方向にモデリングの過程を導きます。

from statsmodels.tsa.deterministic import DeterministicTerm


class BrokenTimeTrend(DeterministicTerm):
    def __init__(self, break_period: int):
        self._break_period = break_period

    def __str__(self):
        return "Broken Time Trend"

    def _eq_attr(self):
        return (self._break_period,)

    def in_sample(self, index: pd.Index):
        nobs = index.shape[0]
        terms = np.zeros((nobs, 2))
        terms[self._break_period :, 0] = 1
        terms[self._break_period :, 1] = np.arange(
            self._break_period + 1, nobs + 1
        )
        return pd.DataFrame(
            terms, columns=["const_break", "trend_break"], index=index
        )

    def out_of_sample(
        self, steps: int, index: pd.Index, forecast_index: pd.Index = None
    ):
        # Always call extend index first
        fcast_index = self._extend_index(index, steps, forecast_index)
        nobs = index.shape[0]
        terms = np.zeros((steps, 2))
        # Assume break period is in-sample
        terms[:, 0] = 1
        terms[:, 1] = np.arange(nobs + 1, nobs + steps + 1)
        return pd.DataFrame(
            terms, columns=["const_break", "trend_break"], index=fcast_index
        )

btt = BrokenTimeTrend(60)
tt = TimeTrend(constant=True, order=1)
index = pd.RangeIndex(100)
det_proc = DeterministicProcess(index, additional_terms=[tt, btt])
det_proc.range(55, 65)
#
const   trend   const_break trend_break
55  1.0 56.0    0.0 0.0
56  1.0 57.0    0.0 0.0
57  1.0 58.0    0.0 0.0
58  1.0 59.0    0.0 0.0
59  1.0 60.0    0.0 0.0
60  1.0 61.0    1.0 61.0
61  1.0 62.0    1.0 62.0
62  1.0 63.0    1.0 63.0
63  1.0 64.0    1.0 64.0
64  1.0 65.0    1.0 65.0
65  1.0 66.0    1.0 66.0

次に、実際の外生的データのシンプルな「ラッパー」を書くことで、サンプル外の外生データを簡単に構築し、予測のために利用します。

class ExogenousProcess(DeterministicTerm):
    def __init__(self, data):
        self._data = data

    def __str__(self):
        return "Custom Exog Process"

    def _eq_attr(self):
        return (id(self._data),)

    def in_sample(self, index: pd.Index):
        return self._data.loc[index]

    def out_of_sample(
        self, steps: int, index: pd.Index, forecast_index: pd.Index = None
    ):
        forecast_index = self._extend_index(index, steps, forecast_index)
        return self._data.loc[forecast_index]

import numpy as np

gen = np.random.default_rng(98765432101234567890)
exog = pd.DataFrame(
    gen.integers(100, size=(300, 2)), columns=["exog1", "exog2"]
)
exog.head()
#
exog1   exog2
0   6   99
1   64  28
2   15  81
3   54  8
4   12  8

DeterministicProcessをサポートするモデル

DeterministicProcessを直接サポートしているモデルはAutoRegだけです。カスタムな項は、決定論的キーワード引数を使用して設定することができます。

注意: カスタムな項を使用するには、 trend="n "とseasonal=Falseが必要なので、すべての決定論的要素はカスタム決定論的項から来ていなければなりません。

いくつかのデータのシミュレーション

ここでは、フーリエ級数によって捕捉された週ごとの季節性を持つデータをシミュレーションしています。

gen = np.random.default_rng(98765432101234567890)
idx = pd.RangeIndex(200)
det_proc = DeterministicProcess(idx, constant=True, period=52, fourier=2)
det_terms = det_proc.in_sample().to_numpy()
params = np.array([1.0, 3, -1, 4, -2])
exog = det_terms @ params
y = np.empty(200)
y[0] = det_terms[0] @ params + gen.standard_normal()
for i in range(1, 200):
    y[i] = 0.9 * y[i - 1] + det_terms[i] @ params + gen.standard_normal()
y = pd.Series(y, index=idx)
ax = y.plot()

image.png

seasonalのデフォルトはFalseですが、deterministicのデフォルトは "c "なので、これを変更する必要があります。

from statsmodels.tsa.api import AutoReg

mod = AutoReg(y, 1, trend="n", deterministic=det_proc)
res = mod.fit()
print(res.summary())
#
AutoReg Model Results                             
==============================================================================
Dep. Variable:                      y   No. Observations:                  200
Model:                     AutoReg(1)   Log Likelihood                -270.964
Method:               Conditional MLE   S.D. of innovations              0.944
Date:                Tue, 23 Feb 2021   AIC                             -0.044
Time:                        13:06:33   BIC                              0.072
Sample:                             1   HQIC                             0.003
                                  200                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.8436      0.172      4.916      0.000       0.507       1.180
sin(1,52)      2.9738      0.160     18.587      0.000       2.660       3.287
cos(1,52)     -0.6771      0.284     -2.380      0.017      -1.235      -0.120
sin(2,52)      3.9951      0.099     40.336      0.000       3.801       4.189
cos(2,52)     -1.7206      0.264     -6.519      0.000      -2.238      -1.203
y.L1           0.9116      0.014     63.264      0.000       0.883       0.940
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            1.0970           +0.0000j            1.0970            0.0000
-----------------------------------------------------------------------------

plot_predictを使用して、予測値とその予測間隔を表示することができます。サンプル外の決定論的な値は、AutoRegに渡された決定論的な処理によって自動的に生成されます。

fig = res.plot_predict(200, 200 + 2 * 52, True)

image.png

auto_reg_forecast = res.predict(200, 211)
auto_reg_forecast
#
200    -3.253482
201    -8.555660
202   -13.607557
203   -18.152622
204   -21.950370
205   -24.790116
206   -26.503171
207   -26.972781
208   -26.141244
209   -24.013773
210   -20.658891
211   -16.205310
dtype: float64
from statsmodels.tsa.api import SARIMAX

det_proc = DeterministicProcess(idx, period=52, fourier=2)
det_terms = det_proc.in_sample()
#
200    -3.253482
201    -8.555660
202   -13.607557
203   -18.152622
204   -21.950370
205   -24.790116
206   -26.503171
207   -26.972781
208   -26.141244
209   -24.013773
210   -20.658891
211   -16.205310
dtype: float64

他のモデルとの併用

他のモデルは、DeterministicProcessを直接サポートしていません。代わりに、外生変数をサポートするモデルに、任意の決定論的項をexogとして手渡すことができます。

外生変数を含むSARIMAXは、SARIMAエラーを含むOLSであることに注意してください。
image.png
決定論的項のパラメータは、式に従って進化する AutoReg と直接比較することはできません。
image.png
xt が決定論的項のみを含む場合、これら 2 つの表現は等価である (MA が存在しないように θ(L)=0 と仮定する)。

mod = SARIMAX(y, order=(1, 0, 0), trend="c", exog=det_terms)
res = mod.fit(disp=False)
print(res.summary())
#
 SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  200
Model:               SARIMAX(1, 0, 0)   Log Likelihood                -293.381
Date:                Tue, 23 Feb 2021   AIC                            600.763
Time:                        13:07:33   BIC                            623.851
Sample:                             0   HQIC                           610.106
                                - 200                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      0.0796      0.140      0.567      0.570      -0.196       0.355
sin(1,52)      9.1916      0.876     10.492      0.000       7.475      10.909
cos(1,52)    -17.4348      0.891    -19.576      0.000     -19.180     -15.689
sin(2,52)      1.2512      0.466      2.683      0.007       0.337       2.165
cos(2,52)    -17.1863      0.434    -39.583      0.000     -18.037     -16.335
ar.L1          0.9957      0.007    150.762      0.000       0.983       1.009
sigma2         1.0748      0.119      9.068      0.000       0.842       1.307
===================================================================================
Ljung-Box (L1) (Q):                   2.16   Jarque-Bera (JB):                 1.03
Prob(Q):                              0.14   Prob(JB):                         0.60
Heteroskedasticity (H):               0.71   Skew:                            -0.14
Prob(H) (two-sided):                  0.16   Kurtosis:                         2.78
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

予測は似ていますが、SARIMAXのパラメータがMLEを用いて推定されているのに対し、AutoRegはOLSを用いて推定されている点が異なります。

sarimax_forecast = res.forecast(12, exog=det_proc.out_of_sample(12))
df = pd.concat([auto_reg_forecast, sarimax_forecast], axis=1)
df.columns = columns = ["AutoReg", "SARIMAX"]
df
#
AutoReg SARIMAX
200 -3.253482   -2.956562
201 -8.555660   -7.985589
202 -13.607557  -12.794073
203 -18.152622  -17.130964
204 -21.950370  -20.760474
205 -24.790116  -23.475512
206 -26.503171  -25.109630
207 -26.972781  -25.546790
208 -26.141244  -24.728386
209 -24.013773  -22.657095
210 -20.658891  -19.397353
211 -16.205310  -15.072386
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
16
Help us understand the problem. What are the problem?