予測、データセットの更新、および「ニュース」¶
このノートブックでは、Statsmodelsを使用して、更新または改訂されたデータセットが標本外予測や欠損データの標本内推定に与える影響を計算する方法について説明します。「Nowcasting」文献のアプローチに従い（参考文献は最後に記載）、状態空間モデルを使用して、データの「ニュース」や影響を計算します。
注意: このノートブックは、Statsmodels v0.12+に適用されます。また、適用されるのは、以下の状態空間モデルまたは関連クラスのみです：sm.tsa.statespace.ExponentialSmoothing、sm.tsa.arima.ARIMA、sm.tsa.SARIMAX、sm.tsa.UnobservedComponents、sm.tsa.VARMAX、およびsm.tsa.DynamicFactor。
[1]:
%matplotlib inline
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
macrodata = sm.datasets.macrodata.load_pandas().data
macrodata.index = pd.period_range('1959Q1', '2009Q3', freq='Q')
予測演習は、多くの場合、モデル選択とパラメータ推定に使用される固定された履歴データのセットから始まります。その後、フィットした選択済みモデル（または複数のモデル）を使用して、標本外予測を作成します。しかし、ほとんどの場合、それで終わりではありません。新しいデータが追加されると、予測誤差を評価し、必要に応じてモデルを更新し、最新の標本外予測を作成する必要があります。このプロセスは「リアルタイム」予測演習と呼ばれることがあります（対照的に、擬似リアルタイム演習は、この手順をシミュレートするものです）。
もし重要なのが予測誤差（MSEなど）に基づく損失関数を最小化することだけであれば、新しいデータが追加された際に、モデル選択、パラメータ推定、標本外予測を完全にやり直すことが望ましい場合があります。この方法を取った場合、新しい予測が変化する理由は2つあります：
新しいデータが得られ、それが新しい情報を提供したこと
予測モデルまたは推定されたパラメータが異なること
このノートブックでは、最初の要因を分離する方法に焦点を当てます。このアプローチは、いわゆる「ナウキャスティング」文献に由来しており、特にBańbura, Giannone, and Reichlin (2011)、Bańbura and Modugno (2014)、およびBańbura et al. (2014)に基づいています。彼らはこの演習を「ニュース」を計算することと説明しており、Statsmodelsでもこの言語を採用しています。
これらの方法は特に多変量モデルで役立つ場合が多いです。なぜなら、複数の変数が同時に更新される可能性があり、どの変数の更新がどの予測変更を引き起こしたのかをすぐに判断することが難しいためです。しかし、単変量モデルの予測修正を考える際にも役立つことがあります。そこで、まずは単純な単変量のケースから説明し、その後で多変量のケースに進むことにします。
改版に関する注意: 私たちが使用しているフレームワークは、新たに観測されたデータポイントによる予測の変化を分解するように設計されています。また、以前に公表されたデータポイントの改版も考慮することができますが、それらを個別に分解することはありません。代わりに、「改版」の総合的な影響のみを示します。
``exog``データに関する注意: 使用しているフレームワークでは、モデル化された変数に対してのみ、新たに観測されたデータポイントによる予測の変化を分解します。これらは、Statsmodelsでendog引数として指定される「左辺」の変数です。このフレームワークでは、exog引数に含まれるようなモデル化されていない「右辺」の変数に対する変化を分解したり考慮したりすることはありません。
単純な単変量の例: AR(1)¶
まず、単純な自己回帰モデルであるAR(1)から始めます：
パラメータ\(\phi\)は時系列の持続性を表します。
このモデルを用いてインフレ率を予測します。
このノートブックで予測更新を簡単に説明するために、平均値が差し引かれたインフレデータを使用します。ただし、実際にはモデルに平均項を追加することは容易です。
[2]:
# De-mean the inflation series
y = macrodata['infl'] - macrodata['infl'].mean()
ステップ 1: 利用可能なデータセットにモデルをフィットさせる¶
ここでは、標本外のシミュレーションを行います。すべてのデータのうち、最後の5つの観測値を除外してモデルを構築およびフィットさせます。これらの値はまだ観測されていないと仮定し、その後のステップで再び分析に追加します。
[3]:
y_pre = y.iloc[:-5]
y_pre.plot(figsize=(15, 3), title='Inflation');
予測を構築するには、まずモデルのパラメータを推定します。これにより、予測を生成するために使用できる結果オブジェクトが返されます。
[4]:
mod_pre = sm.tsa.arima.ARIMA(y_pre, order=(1, 0, 0), trend='n')
res_pre = mod_pre.fit()
print(res_pre.summary())
SARIMAX Results
==============================================================================
Dep. Variable: infl No. Observations: 198
Model: ARIMA(1, 0, 0) Log Likelihood -446.407
Date: Thu, 28 Nov 2024 AIC 896.813
Time: 23:11:44 BIC 903.390
Sample: 03-31-1959 HQIC 899.475
- 06-30-2008
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ar.L1 0.6751 0.043 15.858 0.000 0.592 0.759
sigma2 5.3027 0.367 14.459 0.000 4.584 6.022
===================================================================================
Ljung-Box (L1) (Q): 15.65 Jarque-Bera (JB): 43.04
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 0.85 Skew: 0.18
Prob(H) (two-sided): 0.50 Kurtosis: 5.26
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
結果オブジェクトresから予測を作成するのは簡単です。構築したい予測の数を指定して、forecastメソッドを呼び出すだけです。この場合、標本外の予測を4つ作成します。
[5]:
# Compute the forecasts
forecasts_pre = res_pre.forecast(4)
# Plot the last 3 years of data and the four out-of-sample forecasts
y_pre.iloc[-12:].plot(figsize=(15, 3), label='Data', legend=True)
forecasts_pre.plot(label='Forecast', legend=True);
AR(1)モデルでは、予測を手動で構築することも簡単です。最後に観測された変数を\(y_T\)、\(h\)ステップ先の予測を\(y_{T+h|T}\)とすると、以下のようになります：
ここで、\(\hat \phi\)はAR(1)係数の推定値です。上記のサマリー出力から、この値はモデルの最初のパラメータであり、結果オブジェクトのparams属性から取得できることが分かります。
[6]:
# Get the estimated AR(1) coefficient
phi_hat = res_pre.params[0]
# Get the last observed value of the variable
y_T = y_pre.iloc[-1]
# Directly compute the forecasts at the horizons h=1,2,3,4
manual_forecasts = pd.Series([phi_hat * y_T, phi_hat**2 * y_T,
phi_hat**3 * y_T, phi_hat**4 * y_T],
index=forecasts_pre.index)
# We'll print the two to double-check that they're the same
print(pd.concat([forecasts_pre, manual_forecasts], axis=1))
predicted_mean 0
2008Q3 3.084388 3.084388
2008Q4 2.082323 2.082323
2009Q1 1.405812 1.405812
2009Q2 0.949088 0.949088
C:\Users\user\AppData\Local\Temp\ipykernel_21248\554023716.py:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`
phi_hat = res_pre.params[0]
ステップ2: 新しい観測値から「ニュース」を計算する¶
時間が経過し、新たな観測値が得られたとします。この結果、データセットが拡大し、予測誤差を評価し、次の四半期に向けた更新された予測を作成することができます。
[7]:
# Get the next observation after the "pre" dataset
y_update = y.iloc[-5:-4]
# Print the forecast error
print('Forecast error: %.2f' % (y_update.iloc[0] - forecasts_pre.iloc[0]))
Forecast error: -10.21
更新されたデータセットに基づいて予測を計算するために、appendメソッドを使用して新しい観測データを前のデータセットに追加し、更新された結果オブジェクトres_postを作成します。
デフォルトでは、appendメソッドはモデルのパラメータを再推定しないことに注意してください。これは、ここで求めている通り、新しい情報のみが予測に与える影響を分離したいからです。
[8]:
# Create a new results object by passing the new observations to the `append` method
res_post = res_pre.append(y_update)
# Since we now know the value for 2008Q3, we will only use `res_post` to
# produce forecasts for 2008Q4 through 2009Q2
forecasts_post = pd.concat([y_update, res_post.forecast('2009Q2')])
print(forecasts_post)
2008Q3 -7.121330
2008Q4 -4.807732
2009Q1 -3.245783
2009Q2 -2.191284
Freq: Q-DEC, dtype: float64
この場合、予測誤差は非常に大きく、インフレ率はAR(1)モデルの予測よりも10パーセントポイント以上低かった。（これは主に、世界金融危機前後の原油価格の大きな変動によるものです。）
これをさらに深く分析するために、Statsmodelsを使用して、新しい情報（つまり「ニュース」）が予測に与える影響を分離できます。これは、まだモデルを変更したり、パラメータを再推定したりするつもりはないことを意味します。代わりに、状態空間モデルの結果オブジェクトに利用可能なnewsメソッドを使用します。
Statsmodelsでニュースを計算するには、必ず前の結果オブジェクトまたはデータセットと、更新された結果オブジェクトまたはデータセットが必要です。ここでは、元の結果オブジェクトres_preを前の結果として使用し、先ほど作成したres_postの結果オブジェクトを更新された結果として使用します。
以前の結果と更新された結果のオブジェクトまたはデータセットが揃ったら、news メソッドを呼び出すことで新しい情報を計算できます。ここでは、res_pre.news を呼び出し、最初の引数に更新された結果 res_post を渡します（ただし、2つの結果オブジェクトがある場合、news メソッドはどちらのオブジェクトにも呼び出すことができます）。
比較対象のオブジェクトやデータセットを最初の引数として指定することに加えて、さまざまなその他の引数も受け付けられます。最も重要なのは、考慮したい「影響期間」を指定することです。これらの「影響期間」は、予測された関心のある期間に対応します。つまり、これらの日付は、予測の修正が分解される期間を指定します。
影響期間を指定するには、start、end、および periods の2つを渡す必要があります（Pandasの date_range メソッドに似ています）。もしあなたの時系列がPandasのオブジェクトで、日付または期間のインデックスが関連付けられている場合は、start と end に日付を値として渡すことができます。
[9]:
# Compute the impact of the news on the four periods that we previously
# forecasted: 2008Q3 through 2009Q2
news = res_pre.news(res_post, start='2008Q3', end='2009Q2')
# Note: one alternative way to specify these impact dates is
# `start='2008Q3', periods=4`
変数newsはNewsResultsクラスのオブジェクトであり、res_preと比較してres_postのデータに対する更新の詳細、新しいデータセットにおける新しい情報、およびその新しい情報がstartとendの期間中の予測に与えた影響を含んでいます。
結果を要約する簡単な方法の1つは、summaryメソッドを使用することです。
[10]:
print(news.summary())
News
===============================================================================
Model: ARIMA Original sample: 1959Q1
Date: Thu, 28 Nov 2024 - 2008Q2
Time: 23:11:45 Update through: 2008Q3
# of revisions: 0
# of new datapoints: 1
Impacts for [impacted variable = infl]
=========================================================
impact date estimate (prev) impact of news estimate (new)
---------------------------------------------------------
2008Q3 3.08 -10.21 -7.12
2008Q4 2.08 -6.89 -4.81
2009Q1 1.41 -4.65 -3.25
2009Q2 0.95 -3.14 -2.19
News from updated observations:
===================================================================
update date updated variable observed forecast (prev) news
-------------------------------------------------------------------
2008Q3 infl -7.12 3.08 -10.21
===================================================================
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:1328: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 -7.12
Name: observed, dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
data.iloc[:, 2:] = data.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:1328: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 3.08
Name: forecast (prev), dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
data.iloc[:, 2:] = data.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:1328: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 -10.21
Name: news, dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
data.iloc[:, 2:] = data.iloc[:, 2:].map(
概要出力：このニュース結果オブジェクトのデフォルトのサマリには、次の4つのテーブルが表示されます：
モデルとデータセットのサマリ
更新データからのニュースの詳細
start='2008Q3'からend='2009Q2'までの予測に対する新しい情報の影響のサマリstart='2008Q3'からend='2009Q2'までの予測に対する影響をもたらした更新データの詳細
これらは以下でさらに詳細に説明します。
注記：
この出力を制御するために、
summaryメソッドに渡すことができるいくつかの引数があります。詳細については、ドキュメント / ドキュメンテーション文字列 を確認してください。更新と影響の詳細を示すテーブル（4）は、モデルが多変量であったり、更新が複数回行われたり、影響の日付が多数選ばれたりするとかなり大きくなる可能性があります。デフォルトでは、単変量モデルに対してのみ表示されます。
最初の表：モデルとデータセットの概要
上記の最初の表には次の情報が示されています：
予測が行われたモデルの種類。ここではAR(1)がARIMA(p,d,q)モデルの特別なケースであるため、ARIMAモデルです。
分析が計算された日時。
元の標本期間。ここでは
y_preに対応します。更新された標本期間の終了点。ここでは
y_postの最後の日付です。
第二の表：更新データからのニュース
この表は、更新された標本のなかで更新された観測値に対する前回の結果からの予測を示しています。
注記:
更新されたデータセット
y_postには、以前の観測値に対する修正は含まれていません。もし修正があった場合、その修正ごとの以前の値と更新後の値を示す追加の表が表示されているはずです。
第三の表: 新しい情報の影響の概要
列:
上記の第三の表は、以下を示しています：
各影響日付に対する前回の予測を、「estimate (prev)」列に示しています。
各影響日付に対する新しい情報（「news」）が予測に与えた影響を、「impact of news」列に示しています。
各影響日付に対する更新された予測を、「estimate (new)」列に示しています。
注意事項:
多変量モデルの場合、この表には、各行に関連する影響を受けた変数を説明する追加の列が含まれています。
更新されたデータセット
y_postには、以前観測されたデータポイントに対する修正は含まれていません。もし修正があった場合、その修正が影響日付の予測に与えた影響を示す追加の列がこの表に含まれることになります。estimate (new) = estimate (prev) + impact of newsであることに注意してください。この表は、
summary_impactsメソッドを使用して独立してアクセスすることができます。
私たちの例では:
私たちの例では、表が先に計算した値を示していることに注目してください：
「estimate (prev)」列は、以前のモデルからの予測と一致し、
forecasts_pre変数に含まれています。「estimate (new)」列は、2008年Q3の観測値と2008年Q4〜2009年Q2の更新されたモデルからの予測を含む
forecasts_post変数と一致します。
第4表：更新の詳細とその影響
上記の第4表は、各新しい観測値が各影響日でどのように具体的な影響に変換されたかを示しています。
列の説明：
最初の3列は、関連する更新を記述しています（「更新」は新しい観測値です）：
最初の列（「更新日」）は、更新された変数の日付を示します。
2番目の列（「予測（前回）」）は、前回の結果/データセットに基づいて、更新日で更新変数に対して予測されていた値を示します。
3番目の列（「観測値」）は、更新結果/データセットでその更新変数の実際の観測値を示します。
最後の4列は、特定の更新の影響を記述しています（影響とは「影響期間内での予測の変更」）：
4番目の列（「影響日」）は、指定された更新が影響を与えた日付を示します。
5番目の列（「ニュース」）は、指定された更新に関連する「ニュース」を示します（これは各更新の影響で同じですが、デフォルトでは省略されます）。
6番目の列（「重み」）は、指定された更新の「ニュース」が、影響を受けた変数に与える重みを示します。一般的に、重みは各「更新された変数」/「更新日」/「影響を受けた変数」/「影響日」の組み合わせで異なります。
7番目の列（「影響」）は、指定された更新が、指定された「影響を受けた変数」/「影響日」に与えた影響を示します。
注釈：
多変量モデルでは、この表には各行に対して更新された変数と影響を受けた変数を示す追加の列があります。ここでは1つの変数（「infl」）のみを使用しているため、これらの列は省略されてスペースを節約しています。
デフォルトでは、この表の更新は「スパース化」されており、各行で「更新日」、「予測（前回）」、「観測値」の同じ値が繰り返されるのを避けるために空白が挿入されています。この動作は
sparsify引数を使ってオーバーライドすることができます。影響 = ニュース * 重みであることに注意してください。この表は
summary_detailsメソッドを使用して独立してアクセスすることができます。
例：
2008年Q3の更新と影響日2008年Q3において、重みは1です。これは、1つの変数しかなく、2008年Q3のデータを組み込んだ時点で、この日付の「予測」に関して残る曖昧さがないためです。したがって、2008年Q3でのこの変数に関するすべての「ニュース」は直接「予測」に反映されます。
補遺: ニュース、重み、影響の手動計算¶
この単純な一変量モデルの例では、上記に示したすべての値を手動で計算するのは簡単です。まず、予測の式 \(y_{T+h|T} = \phi^h y_T\) を思い出し、次に \(y_{T+h|T+1} = \phi^h y_{T+1}\) が成り立つことに注意してください。最後に、\(y_{T|T+1} = y_T\) であることに注意してください。これは、\(T+1\) までの観測値が分かっていれば、\(y_T\) の値も分かるからです。
ニュース: 「ニュース」とは、単に新しい観測に関連する予測誤差です。したがって、観測 \(T+1\) に関連するニュースは次のようになります：
影響: ニュースの影響とは、更新された予測と以前の予測との差であり、\(i_h \equiv y_{T+h|T+1} - y_{T+h|T}\) です。
\(h=1, \dots, 4\) の場合の以前の予測は次の通りです: \(\begin{pmatrix} \phi y_T & \phi^2 y_T & \phi^3 y_T & \phi^4 y_T \end{pmatrix}'\).
\(h=1, \dots, 4\) の場合の更新された予測は次の通りです: \(\begin{pmatrix} y_{T+1} & \phi y_{T+1} & \phi^2 y_{T+1} & \phi^3 y_{T+1} \end{pmatrix}'\).
したがって、影響は次のようになります：
重み: 重みを計算するには、影響を予測誤差 \(n_{T+1}\) の観点で書き換えることができることに注意すれば十分です。
したがって、重みは単に次のようになります：\(w = \begin{pmatrix} 1 \\ \phi \\ \phi^2 \\ \phi^3 \end{pmatrix}\)
これが news メソッドが計算したものであることを確認できます。
[11]:
# Print the news, computed by the `news` method
print(news.news)
# Manually compute the news
print()
print((y_update.iloc[0] - phi_hat * y_pre.iloc[-1]).round(6))
update date updated variable
2008Q3 infl -10.205718
Name: news, dtype: float64
-10.205718
[12]:
# Print the total impacts, computed by the `news` method
# (Note: news.total_impacts = news.revision_impacts + news.update_impacts, but
# here there are no data revisions, so total and update impacts are the same)
print(news.total_impacts)
# Manually compute the impacts
print()
print(forecasts_post - forecasts_pre)
impacted variable infl
impact date
2008Q3 -10.205718
2008Q4 -6.890055
2009Q1 -4.651595
2009Q2 -3.140371
2008Q3 -10.205718
2008Q4 -6.890055
2009Q1 -4.651595
2009Q2 -3.140371
Freq: Q-DEC, dtype: float64
[13]:
# Print the weights, computed by the `news` method
print(news.weights)
# Manually compute the weights
print()
print(np.array([1, phi_hat, phi_hat**2, phi_hat**3]).round(6))
impact date 2008Q3 2008Q4 2009Q1 2009Q2
impacted variable infl infl infl infl
update date updated variable
2008Q3 infl 1.0 0.675117 0.455783 0.307707
[1. 0.675117 0.455783 0.307707]
多変量の例: 動的因子¶
この例では、個人消費支出（PCE）価格指数と消費者物価指数（CPI）を使用して、月次のコア価格インフレを予測するために動的因子モデルを考えます。これらの指標はどちらも米国経済の価格を追跡しており、同様の元データに基づいていますが、定義にいくつかの違いがあります。それにもかかわらず、これらは互いに比較的よく連動しているため、単一の動的因子を用いてそれらを共同でモデル化するのは合理的に思えます。
このアプローチが有用である理由の一つは、CPIがPCEよりも月初に発表されることです。したがって、CPIが発表され次第、その追加データポイントを使って動的因子モデルを更新し、その月のPCE発表に対する予測を改善することができます。この種の分析のより詳細なバージョンは、KnotekとZaman（2017）で利用可能です。
まず、FREDからコアCPIとPCE物価指数データをダウンロードし、それらを年率換算の月次インフレ率に変換し、2つの外れ値を除去し、各系列の平均値を引いて調整します（動的因子モデルでは）。
[14]:
import pandas_datareader as pdr
levels = pdr.get_data_fred(['PCEPILFE', 'CPILFESL'], start='1999', end='2019').to_period('M')
infl = np.log(levels).diff().iloc[1:] * 1200
infl.columns = ['PCE', 'CPI']
# Remove two outliers and de-mean the series
infl['PCE'].loc['2001-09':'2001-10'] = np.nan
C:\Users\user\AppData\Local\Temp\ipykernel_21248\2462809934.py:7: FutureWarning: ChainedAssignmentError: behaviour will change in pandas 3.0!
You are setting values through chained assignment. Currently this works in certain cases, but when using Copy-on-Write (which will become the default behaviour in pandas 3.0) this will never work to update the original DataFrame or Series, because the intermediate object on which we are setting values will behave as a copy.
A typical example is when you are setting values in a column of a DataFrame, like:
df["col"][row_indexer] = value
Use `df.loc[row_indexer, "col"] = values` instead, to perform the assignment in a single step and ensure this keeps updating the original `df`.
See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
infl['PCE'].loc['2001-09':'2001-10'] = np.nan
これがどのように機能するかを示すために、2017年4月14日、つまり2017年3月のCPI発表の日を想定します。複数の更新の影響を一度に示すために、1月末以来データを更新していないと仮定します。したがって：
以前のデータセットは、2017年1月までのPCEとCPIのすべての値で構成されます。
更新されたデータセットは、さらに2017年2月と3月のCPI、および2017年2月のPCEデータを含みます。ただし、2017年3月のPCE（PCE価格指数は2017年5月1日まで発表されなかった）はまだ含まれていません。
[15]:
# Previous dataset runs through 2017-02
y_pre = infl.loc[:'2017-01'].copy()
const_pre = np.ones(len(y_pre))
print(y_pre.tail())
PCE CPI
DATE
2016-09 1.591601 2.022262
2016-10 1.540990 1.445830
2016-11 0.533425 1.631694
2016-12 1.393060 2.109728
2017-01 3.203951 2.623570
[16]:
# For the updated dataset, we'll just add in the
# CPI value for 2017-03
y_post = infl.loc[:'2017-03'].copy()
y_post.loc['2017-03', 'PCE'] = np.nan
const_post = np.ones(len(y_post))
# Notice the missing value for PCE in 2017-03
print(y_post.tail())
PCE CPI
DATE
2016-11 0.533425 1.631694
2016-12 1.393060 2.109728
2017-01 3.203951 2.623570
2017-02 2.123190 2.541355
2017-03 NaN -0.258197
この特定の例を選んだ理由は、2017年3月にコアCPI（消費者物価指数）が2010年以来初めて低下し、この情報がその月のコアPCE（個人消費支出）価格予測に役立つ可能性があるからです。下のグラフは、2017年4月14日に観測されたであろうCPIとPCEの価格データを示しています。
\(\dagger\) この記述は完全に正確ではありません。なぜなら、CPIとPCE価格指数は事後にある程度修正されることがあるためです。その結果、私たちが取得しているシリーズは、2017年4月14日に観測されたものとは正確に一致しません。この問題は、FREDではなく、ALFREDからアーカイブされたデータを取得することで修正できますが、このチュートリアルでは提供されているデータで十分です。
[17]:
# Plot the updated dataset
fig, ax = plt.subplots(figsize=(15, 3))
y_post.plot(ax=ax)
ax.hlines(0, '2009', '2017-06', linewidth=1.0)
ax.set_xlim('2009', '2017-06');
この演習を行うために、まずDynamicFactorモデルを構築してフィットさせます。具体的には：
単一の動的因子を使用しています（
k_factors=1）。因子の動態をAR(6)モデルでモデル化しています（
factor_order=6）。使用しているインフレシリーズが平均ゼロではないため、外生変数として1のベクトル（
exog=const_pre）を含めています。
[18]:
mod_pre = sm.tsa.DynamicFactor(y_pre, exog=const_pre, k_factors=1, factor_order=6)
res_pre = mod_pre.fit()
print(res_pre.summary())
Statespace Model Results
=============================================================================================
Dep. Variable: ['PCE', 'CPI'] No. Observations: 216
Model: DynamicFactor(factors=1, order=6) Log Likelihood -522.884
+ 1 regressors AIC 1069.768
Date: Thu, 28 Nov 2024 BIC 1110.271
Time: 23:18:58 HQIC 1086.131
Sample: 02-28-1999
- 01-31-2017
Covariance Type: opg
===================================================================================
Ljung-Box (L1) (Q): 4.50, 0.54 Jarque-Bera (JB): 13.09, 12.63
Prob(Q): 0.03, 0.46 Prob(JB): 0.00, 0.00
Heteroskedasticity (H): 0.56, 0.44 Skew: 0.18, -0.16
Prob(H) (two-sided): 0.02, 0.00 Kurtosis: 4.15, 4.14
Results for equation PCE
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
loading.f1 0.5499 0.061 9.037 0.000 0.431 0.669
beta.const 1.7041 0.095 17.960 0.000 1.518 1.890
Results for equation CPI
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
loading.f1 0.9033 0.102 8.876 0.000 0.704 1.103
beta.const 1.9623 0.137 14.357 0.000 1.694 2.230
Results for factor equation f1
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
L1.f1 0.1247 0.069 1.806 0.071 -0.011 0.260
L2.f1 0.1823 0.072 2.548 0.011 0.042 0.323
L3.f1 0.0178 0.073 0.244 0.807 -0.125 0.160
L4.f1 -0.0700 0.078 -0.893 0.372 -0.224 0.084
L5.f1 0.1561 0.068 2.308 0.021 0.024 0.289
L6.f1 0.1376 0.075 1.842 0.065 -0.009 0.284
Error covariance matrix
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
sigma2.PCE 0.5422 0.065 8.287 0.000 0.414 0.670
sigma2.CPI 1.15e-08 0.144 7.96e-08 1.000 -0.283 0.283
==============================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
フィットしたモデルを手に入れたので、次に2017年3月のCPI観測に関連するニュースと影響を構築します。更新されたデータは2017年2月と2017年3月の一部で、3月と4月の影響を検討します。
一変量モデルの例では、まず更新された結果オブジェクトを作成し、それを news メソッドに渡しました。ここでは、更新されたデータセットを直接渡してニュースを作成しています。
以下の点に注意してください：
y_postには更新されたデータセット全体が含まれています（新しいデータポイントだけではありません）。更新された
exog配列も渡す必要がありました。この配列は両方をカバーしなければなりません：y_postに関連する期間全体y_postの終了後から最後の影響日までの追加のデータポイント（endで指定）
ここでは、y_post は2017年3月で終了するため、exog はさらに1期間延長して2017年4月までをカバーする必要がありました。
[19]:
# Create the news results
# Note
const_post_plus1 = np.ones(len(y_post) + 1)
news = res_pre.news(y_post, exog=const_post_plus1, start='2017-03', end='2017-04')
注:
上記の一変量の例では、まず新しい結果オブジェクトを作成し、それを
newsメソッドに渡しました。ここでも同様のことができますが、追加の手順が必要です。y_postの終了後の期間について影響を要求しているため、その期間のexog変数の追加値もnewsに渡す必要があります：res_post = res_pre.apply(y_post, exog=const_post) news = res_pre.news(res_post, exog=[1.], start='2017-03', end='2017-04')
news を計算したので、summary を表示することが結果を確認する便利な方法です。
[20]:
# Show the summary of the news results
print(news.summary())
News
===============================================================================
Model: DynamicFactor Original sample: 1999-02
Date: Thu, 28 Nov 2024 - 2017-01
Time: 23:18:58 Update through: 2017-03
# of revisions: 0
# of new datapoints: 3
Impacts
===========================================================================
impact date impacted variable estimate (prev) impact of news estimate (new)
---------------------------------------------------------------------------
2017-03 CPI 2.07 -2.33 -0.26
NaT PCE 1.77 -1.42 0.35
2017-04 CPI 1.90 -0.23 1.67
NaT PCE 1.67 -0.14 1.53
News from updated observations:
===================================================================
update date updated variable observed forecast (prev) news
-------------------------------------------------------------------
2017-02 CPI 2.54 2.24 0.30
PCE 2.12 1.87 0.25
2017-03 CPI -0.26 2.07 -2.33
===================================================================
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:936: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 2.07
1 1.77
2 1.90
3 1.67
Name: estimate (prev), dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
impacts.iloc[:, 2:] = impacts.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:936: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 0.00
1 0.00
2 0.00
3 0.00
Name: impact of revisions, dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
impacts.iloc[:, 2:] = impacts.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:936: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 -2.33
1 -1.42
2 -0.23
3 -0.14
Name: impact of news, dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
impacts.iloc[:, 2:] = impacts.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:936: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 -2.33
1 -1.42
2 -0.23
3 -0.14
Name: total impact, dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
impacts.iloc[:, 2:] = impacts.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:936: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 -0.26
1 0.35
2 1.67
3 1.53
Name: estimate (new), dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
impacts.iloc[:, 2:] = impacts.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:1328: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 2.54
1 2.12
2 -0.26
Name: observed, dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
data.iloc[:, 2:] = data.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:1328: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 2.24
1 1.87
2 2.07
Name: forecast (prev), dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
data.iloc[:, 2:] = data.iloc[:, 2:].map(
C:\Users\user\projects\statsmodels-main\statsmodels\tsa\statespace\news.py:1328: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise in a future error of pandas. Value '0 0.30
1 0.25
2 -2.33
Name: news, dtype: object' has dtype incompatible with float64, please explicitly cast to a compatible dtype first.
data.iloc[:, 2:] = data.iloc[:, 2:].map(
複数の変数があるため、デフォルトではサマリーは更新データからのニュースと総影響のみを表示します。
最初の表から、更新されたデータセットには3つの新しいデータポイントが含まれており、そのほとんどの「ニュース」は2017年3月の非常に低い値から来ていることが分かります。
2番目の表は、これらの3つのデータポイントが2017年3月のPCEの推定値に大きな影響を与えたことを示しています（まだ観測されていない）。この推定値はほぼ1.5パーセントポイント下方修正されました。
更新データは、最初の標本外月である2017年4月の予測にも影響を与えました。新しいデータを組み込んだ後、この月のCPIとPCEインフレの予測はそれぞれ0.29ポイントと0.17ポイント下方修正されました。
これらのテーブルは「ニュース」と総影響を示していますが、各影響がどの更新されたデータポイントによって引き起こされたかは示していません。その情報を見るためには、詳細テーブルを確認する必要があります。
詳細テーブルを見る一つの方法は、summary メソッドに include_details=True を渡すことです。ただし、上記のテーブルを繰り返さないように、直接 summary_details メソッドを呼び出すことにします。
[21]:
print(news.summary_details())
Details of news
======================================================================================================================
update date updated variable observed forecast (prev) impact date impacted variable news weight impact
----------------------------------------------------------------------------------------------------------------------
2017-02 CPI 2.54 2.24 2017-03 CPI 0.30 0.00 0.00
PCE 0.30 0.00 0.00
2017-04 CPI 0.30 0.18 0.06
PCE 0.30 0.11 0.03
PCE 2.12 1.87 CPI 0.25 0.00 0.00
PCE 0.25 0.00 0.00
2017-03 CPI -0.26 2.07 2017-03 CPI -2.33 1.00 -2.33
PCE -2.33 0.61 -1.42
2017-04 CPI -2.33 0.12 -0.29
PCE -2.33 0.08 -0.18
======================================================================================================================
この表は、上記で説明した2017年4月のPCE予測の改訂のほとんどが、2017年3月のCPI発表に関連するニュースから来ていることを示しています。それに対して、2017年2月のCPI発表は4月の予測にほとんど影響を与えず、2月のPCE発表は実質的に影響を与えませんでした。
参考文献¶
Bańbura, Marta, Domenico Giannone, and Lucrezia Reichlin. 「Nowcasting.」 The Oxford Handbook of Economic Forecasting. 2011年7月8日.
Bańbura, Marta, Domenico Giannone, Michele Modugno, and Lucrezia Reichlin. 「Now-casting and the real-time data flow.」 In Handbook of economic forecasting, vol. 2, pp. 195-237. Elsevier, 2013.
Bańbura, Marta, and Michele Modugno. 「Maximum likelihood estimation of factor models on datasets with arbitrary pattern of missing data.」 Journal of Applied Econometrics 29, no. 1 (2014): 133-160.
Knotek, Edward S., and Saeed Zaman. 「Nowcasting US headline and core inflation.」 Journal of Money, Credit and Banking 49, no. 5 (2017): 931-968.