状態空間モデル：局所線形トレンド¶

このノートブックでは、statsmodelsの状態空間クラスを拡張してカスタムモデルを作成し、推定する方法を説明します。ここでは、ローカルリニアトレンドモデルを開発します。

ローカルリニアトレンドモデルは次の形式を持ちます（記号や詳細についてはDurbinとKoopman 2012年の第3章2節を参照）：

\[\begin{split}\begin{align} y_t & = \mu_t + \varepsilon_t \qquad & \varepsilon_t \sim N(0, \sigma_\varepsilon^2) \\ \mu_{t+1} & = \mu_t + \nu_t + \xi_t & \xi_t \sim N(0, \sigma_\xi^2) \\ \nu_{t+1} & = \nu_t + \zeta_t & \zeta_t \sim N(0, \sigma_\zeta^2) \end{align}\end{split}\]

これが状態空間形式に変換できることは簡単に分かります：

\[\begin{split}\begin{align} y_t & = \begin{pmatrix} 1 & 0 \end{pmatrix} \begin{pmatrix} \mu_t \\ \nu_t \end{pmatrix} + \varepsilon_t \\ \begin{pmatrix} \mu_{t+1} \\ \nu_{t+1} \end{pmatrix} & = \begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix} \begin{pmatrix} \mu_t \\ \nu_t \end{pmatrix} + \begin{pmatrix} \xi_t \\ \zeta_t \end{pmatrix} \end{align}\end{split}\]

状態空間表現の多くは既知の値で構成されていることに注意してください。実際、推定すべきパラメータが現れるのは分散／共分散行列の部分だけです：

\[\begin{split}\begin{align} H_t & = \begin{bmatrix} \sigma_\varepsilon^2 \end{bmatrix} \\ Q_t & = \begin{bmatrix} \sigma_\xi^2 & 0 \\ 0 & \sigma_\zeta^2 \end{bmatrix} \end{align}\end{split}\]

[1]:

%matplotlib inline

import numpy as np
import pandas as pd
from scipy.stats import norm
import statsmodels.api as sm
import matplotlib.pyplot as plt

既存のインフラストラクチャ（カルマンフィルタリングや最尤推定を含む）を活用するために、statsmodels.tsa.statespace.MLEModelから拡張する新しいクラスを作成します。以下の項目を指定する必要があります：

k_states、k_posdef: これらの2つのパラメータは、初期化時に基底クラスに提供する必要があります。これらはそれぞれ、状態ベクトル（\(\begin{pmatrix} \mu_t & \nu_t \end{pmatrix}'\)）と状態誤差ベクトル（\(\begin{pmatrix} \xi_t & \zeta_t \end{pmatrix}'\)）のサイズを、状態空間モデルに通知します。内生変数ベクトルの次元は、endog配列から推定できるため、指定する必要はありません。
update: updateメソッド（引数paramsを取る）は指定する必要があります（fit()が呼び出される際に最尤推定を計算するために使用されます）。このメソッドはパラメータを受け取り、それらを適切な状態空間行列に埋め込みます。例えば、以下のように、paramsベクトルは分散パラメータ（\(\begin{pmatrix} \sigma_\varepsilon^2 & \sigma_\xi^2 & \sigma_\zeta^2 \end{pmatrix}\)）を含んでおり、updateメソッドはそれらを観測および状態共分散行列に配置する必要があります。一般的に、パラメータベクトルはすべての状態空間行列のさまざまな場所にマッピングされることがあります。
状態空間行列: デフォルトでは、すべての状態空間行列（obs_intercept, design, obs_cov, state_intercept, transition, selection, state_cov）はゼロに設定されます。固定値（ここでの設計行列や遷移行列のようなもの）は初期化時に設定できますが、パラメータに依存する値はupdateメソッドで設定する必要があります。選択行列（ここでは単位行列）の設定を忘れがちですが、設定しないと遷移方程式に確率的成分がない非常に異なるモデルになります。
開始パラメータ: 開始パラメータは設定する必要があります。たとえゼロベクトルであっても、良い開始パラメータはデータから見つけることができます。勾配法による最尤推定（ここで使用されている方法）は開始パラメータに敏感であるため、可能であれば良いものを選択することが重要です。ここではあまり問題になりませんが（分散としてゼロを設定するべきではありません）。
初期化: 定義された状態空間行列に加えて、すべての状態空間モデルは状態ベクトルの初期分布の平均と分散で初期化する必要があります。分布が既知であれば、initialize_known(initial_state, initial_state_cov)を呼び出すことができ、モデルが定常的であれば（例えばARMAモデルなど）、initialize_stationaryを使用できます。それ以外の場合は、initialize_approximate_diffuseが合理的な汎用初期化となります（正確な拡散初期化はまだ利用できません）。ローカル線形トレンドモデルは定常ではなく（ランダムウォークで構成されているため）、分布が一般に知られていないため、以下でinitialize_approximate_diffuseを使用します。

上記は成功するモデルに必要な最小限の項目です。また、必ずしも設定しなくてもよいが、いくつかのアプリケーションには役立つ、あるいは重要な項目もあります：

transform / untransform: fitが呼び出されると、バックグラウンドで最適化アルゴリズムが勾配法を使用して最尤関数を最大化するパラメータを選択します。デフォルトでは無制約の最適化が使用されるため、任意のパラメータ値を選択することができます。多くの場合、これは望ましくない挙動です。例えば、分散は負であってはならないため、これを回避するために、transformメソッドは最適化アルゴリズムが提供する制約のないパラメータベクトルを受け取り、最尤評価に使用する制約付きパラメータベクトルを返します。untransformはその逆の操作を行います。
param_names: この内部メソッドを使用して、推定されたパラメータに名前を付けることができます。これにより、例えばサマリーが意味のある名前を提供するようになります。名前が指定されていない場合、パラメータはparam0、param1などと名付けられます。

[2]:

"""
Univariate Local Linear Trend Model
"""
class LocalLinearTrend(sm.tsa.statespace.MLEModel):
    def __init__(self, endog):
        # Model order
        k_states = k_posdef = 2

        # Initialize the statespace
        super(LocalLinearTrend, self).__init__(
            endog, k_states=k_states, k_posdef=k_posdef,
            initialization='approximate_diffuse',
            loglikelihood_burn=k_states
        )

        # Initialize the matrices
        self.ssm['design'] = np.array([1, 0])
        self.ssm['transition'] = np.array([[1, 1],
                                       [0, 1]])
        self.ssm['selection'] = np.eye(k_states)

        # Cache some indices
        self._state_cov_idx = ('state_cov',) + np.diag_indices(k_posdef)

    @property
    def param_names(self):
        return ['sigma2.measurement', 'sigma2.level', 'sigma2.trend']

    @property
    def start_params(self):
        return [np.std(self.endog)]*3

    def transform_params(self, unconstrained):
        return unconstrained**2

    def untransform_params(self, constrained):
        return constrained**0.5

    def update(self, params, *args, **kwargs):
        params = super(LocalLinearTrend, self).update(params, *args, **kwargs)

        # Observation covariance
        self.ssm['obs_cov',0,0] = params[0]

        # State covariance
        self.ssm[self._state_cov_idx] = params[1:]

この簡単なモデルを使用することで、ローカル線形トレンドモデルからパラメータを推定することができます。以下の例は、Commandeur and Koopman (2007) のセクション3.4からのもので、フィンランドにおける自動車死亡事故のモデル化です。

[3]:

import requests
from io import BytesIO
from zipfile import ZipFile

# Download the dataset
ck = requests.get('http://staff.feweb.vu.nl/koopman/projects/ckbook/OxCodeAll.zip').content
zipped = ZipFile(BytesIO(ck))
df = pd.read_table(
    BytesIO(zipped.read('OxCodeIntroStateSpaceBook/Chapter_2/NorwayFinland.txt')),
    skiprows=1, header=None, sep='\s+', engine='python',
    names=['date','nf', 'ff']
)

---------------------------------------------------------------------------
gaierror                                  Traceback (most recent call last)
File ~\projects\statsmodels\venv\lib\site-packages\urllib3\connection.py:196, in HTTPConnection._new_conn(self)
    195 try:
--> 196     sock = connection.create_connection(
    197         (self._dns_host, self.port),
    198         self.timeout,
    199         source_address=self.source_address,
    200         socket_options=self.socket_options,
    201     )
    202 except socket.gaierror as e:

File ~\projects\statsmodels\venv\lib\site-packages\urllib3\util\connection.py:60, in create_connection(address, timeout, source_address, socket_options)
     58     raise LocationParseError(f"'{host}', label empty or too long") from None
---> 60 for res in socket.getaddrinfo(host, port, family, socket.SOCK_STREAM):
     61     af, socktype, proto, canonname, sa = res

File c:\users\user\appdata\local\programs\python\python310\lib\socket.py:955, in getaddrinfo(host, port, family, type, proto, flags)
    954 addrlist = []
--> 955 for res in _socket.getaddrinfo(host, port, family, type, proto, flags):
    956     af, socktype, proto, canonname, sa = res

gaierror: [Errno 11001] getaddrinfo failed

The above exception was the direct cause of the following exception:

NameResolutionError                       Traceback (most recent call last)
File ~\projects\statsmodels\venv\lib\site-packages\urllib3\connectionpool.py:789, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
    788 # Make the request on the HTTPConnection object
--> 789 response = self._make_request(
    790     conn,
    791     method,
    792     url,
    793     timeout=timeout_obj,
    794     body=body,
    795     headers=headers,
    796     chunked=chunked,
    797     retries=retries,
    798     response_conn=response_conn,
    799     preload_content=preload_content,
    800     decode_content=decode_content,
    801     **response_kw,
    802 )
    804 # Everything went great!

File ~\projects\statsmodels\venv\lib\site-packages\urllib3\connectionpool.py:495, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
    494 try:
--> 495     conn.request(
    496         method,
    497         url,
    498         body=body,
    499         headers=headers,
    500         chunked=chunked,
    501         preload_content=preload_content,
    502         decode_content=decode_content,
    503         enforce_content_length=enforce_content_length,
    504     )
    506 # We are swallowing BrokenPipeError (errno.EPIPE) since the server is
    507 # legitimately able to close the connection after sending a valid response.
    508 # With this behaviour, the received response is still readable.

File ~\projects\statsmodels\venv\lib\site-packages\urllib3\connection.py:398, in HTTPConnection.request(self, method, url, body, headers, chunked, preload_content, decode_content, enforce_content_length)
    397     self.putheader(header, value)
--> 398 self.endheaders()
    400 # If we're given a body we start sending that in chunks.

File c:\users\user\appdata\local\programs\python\python310\lib\http\client.py:1277, in HTTPConnection.endheaders(self, message_body, encode_chunked)
   1276     raise CannotSendHeader()
-> 1277 self._send_output(message_body, encode_chunked=encode_chunked)

File c:\users\user\appdata\local\programs\python\python310\lib\http\client.py:1037, in HTTPConnection._send_output(self, message_body, encode_chunked)
   1036 del self._buffer[:]
-> 1037 self.send(msg)
   1039 if message_body is not None:
   1040
   1041     # create a consistent interface to message_body

File c:\users\user\appdata\local\programs\python\python310\lib\http\client.py:975, in HTTPConnection.send(self, data)
    974 if self.auto_open:
--> 975     self.connect()
    976 else:

File ~\projects\statsmodels\venv\lib\site-packages\urllib3\connection.py:236, in HTTPConnection.connect(self)
    235 def connect(self) -> None:
--> 236     self.sock = self._new_conn()
    237     if self._tunnel_host:
    238         # If we're tunneling it means we're connected to our proxy.

File ~\projects\statsmodels\venv\lib\site-packages\urllib3\connection.py:203, in HTTPConnection._new_conn(self)
    202 except socket.gaierror as e:
--> 203     raise NameResolutionError(self.host, self, e) from e
    204 except SocketTimeout as e:

NameResolutionError: <urllib3.connection.HTTPConnection object at 0x000001FA720C1A20>: Failed to resolve 'personal.vu.nl' ([Errno 11001] getaddrinfo failed)

The above exception was the direct cause of the following exception:

MaxRetryError                             Traceback (most recent call last)
File ~\projects\statsmodels\venv\lib\site-packages\requests\adapters.py:667, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
    666 try:
--> 667     resp = conn.urlopen(
    668         method=request.method,
    669         url=url,
    670         body=request.body,
    671         headers=request.headers,
    672         redirect=False,
    673         assert_same_host=False,
    674         preload_content=False,
    675         decode_content=False,
    676         retries=self.max_retries,
    677         timeout=timeout,
    678         chunked=chunked,
    679     )
    681 except (ProtocolError, OSError) as err:

File ~\projects\statsmodels\venv\lib\site-packages\urllib3\connectionpool.py:843, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
    841     new_e = ProtocolError("Connection aborted.", new_e)
--> 843 retries = retries.increment(
    844     method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
    845 )
    846 retries.sleep()

File ~\projects\statsmodels\venv\lib\site-packages\urllib3\util\retry.py:519, in Retry.increment(self, method, url, response, error, _pool, _stacktrace)
    518     reason = error or ResponseError(cause)
--> 519     raise MaxRetryError(_pool, url, reason) from reason  # type: ignore[arg-type]
    521 log.debug("Incremented Retry for (url='%s'): %r", url, new_retry)

MaxRetryError: HTTPConnectionPool(host='personal.vu.nl', port=80): Max retries exceeded with url: /s.j.koopman/projects/ckbook/OxCodeAll.zip (Caused by NameResolutionError("<urllib3.connection.HTTPConnection object at 0x000001FA720C1A20>: Failed to resolve 'personal.vu.nl' ([Errno 11001] getaddrinfo failed)"))

During handling of the above exception, another exception occurred:

ConnectionError                           Traceback (most recent call last)
Cell In[3], line 6
      3 from zipfile import ZipFile
      5 # Download the dataset
----> 6 ck = requests.get('http://staff.feweb.vu.nl/koopman/projects/ckbook/OxCodeAll.zip').content
      7 zipped = ZipFile(BytesIO(ck))
      8 df = pd.read_table(
      9     BytesIO(zipped.read('OxCodeIntroStateSpaceBook/Chapter_2/NorwayFinland.txt')),
     10     skiprows=1, header=None, sep='\s+', engine='python',
     11     names=['date','nf', 'ff']
     12 )

File ~\projects\statsmodels\venv\lib\site-packages\requests\api.py:73, in get(url, params, **kwargs)
     62 def get(url, params=None, **kwargs):
     63     r"""Sends a GET request.
     64
     65     :param url: URL for the new :class:`Request` object.
   (...)
     70     :rtype: requests.Response
     71     """
---> 73     return request("get", url, params=params, **kwargs)

File ~\projects\statsmodels\venv\lib\site-packages\requests\api.py:59, in request(method, url, **kwargs)
     55 # By using the 'with' statement we are sure the session is closed, thus we
     56 # avoid leaving sockets open which can trigger a ResourceWarning in some
     57 # cases, and look like a memory leak in others.
     58 with sessions.Session() as session:
---> 59     return session.request(method=method, url=url, **kwargs)

File ~\projects\statsmodels\venv\lib\site-packages\requests\sessions.py:589, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
    584 send_kwargs = {
    585     "timeout": timeout,
    586     "allow_redirects": allow_redirects,
    587 }
    588 send_kwargs.update(settings)
--> 589 resp = self.send(prep, **send_kwargs)
    591 return resp

File ~\projects\statsmodels\venv\lib\site-packages\requests\sessions.py:724, in Session.send(self, request, **kwargs)
    721 if allow_redirects:
    722     # Redirect resolving generator.
    723     gen = self.resolve_redirects(r, request, **kwargs)
--> 724     history = [resp for resp in gen]
    725 else:
    726     history = []

File ~\projects\statsmodels\venv\lib\site-packages\requests\sessions.py:724, in <listcomp>(.0)
    721 if allow_redirects:
    722     # Redirect resolving generator.
    723     gen = self.resolve_redirects(r, request, **kwargs)
--> 724     history = [resp for resp in gen]
    725 else:
    726     history = []

File ~\projects\statsmodels\venv\lib\site-packages\requests\sessions.py:265, in SessionRedirectMixin.resolve_redirects(self, resp, req, stream, timeout, verify, cert, proxies, yield_requests, **adapter_kwargs)
    263     yield req
    264 else:
--> 265     resp = self.send(
    266         req,
    267         stream=stream,
    268         timeout=timeout,
    269         verify=verify,
    270         cert=cert,
    271         proxies=proxies,
    272         allow_redirects=False,
    273         **adapter_kwargs,
    274     )
    276     extract_cookies_to_jar(self.cookies, prepared_request, resp.raw)
    278     # extract redirect url, if any, for the next loop

File ~\projects\statsmodels\venv\lib\site-packages\requests\sessions.py:703, in Session.send(self, request, **kwargs)
    700 start = preferred_clock()
    702 # Send the request
--> 703 r = adapter.send(request, **kwargs)
    705 # Total elapsed time of the request (approximately)
    706 elapsed = preferred_clock() - start

File ~\projects\statsmodels\venv\lib\site-packages\requests\adapters.py:700, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
    696     if isinstance(e.reason, _SSLError):
    697         # This branch is for urllib3 v1.22 and later.
    698         raise SSLError(e, request=request)
--> 700     raise ConnectionError(e, request=request)
    702 except ClosedPoolError as e:
    703     raise ConnectionError(e, request=request)

ConnectionError: HTTPConnectionPool(host='personal.vu.nl', port=80): Max retries exceeded with url: /s.j.koopman/projects/ckbook/OxCodeAll.zip (Caused by NameResolutionError("<urllib3.connection.HTTPConnection object at 0x000001FA720C1A20>: Failed to resolve 'personal.vu.nl' ([Errno 11001] getaddrinfo failed)"))

Since we defined the local linear trend model as extending from MLEModel, the fit() method is immediately available, just as in other statsmodels maximum likelihood classes. Similarly, the returned results class supports many of the same post-estimation results, like the summary method.

[ ]:

# Load Dataset
df.index = pd.date_range(start='%d-01-01' % df.date[0], end='%d-01-01' % df.iloc[-1, 0], freq='AS')

# Log transform
df['lff'] = np.log(df['ff'])

# Setup the model
mod = LocalLinearTrend(df['lff'])

# Fit it using MLE (recall that we are fitting the three variance parameters)
res = mod.fit(disp=False)
print(res.summary())

Finally, we can do post-estimation prediction and forecasting. Notice that the end period can be specified as a date.

[ ]:

# Perform prediction and forecasting
predict = res.get_prediction()
forecast = res.get_forecast('2014')

[ ]:

fig, ax = plt.subplots(figsize=(10,4))

# Plot the results
df['lff'].plot(ax=ax, style='k.', label='Observations')
predict.predicted_mean.plot(ax=ax, label='One-step-ahead Prediction')
predict_ci = predict.conf_int(alpha=0.05)
predict_index = np.arange(len(predict_ci))
ax.fill_between(predict_index[2:], predict_ci.iloc[2:, 0], predict_ci.iloc[2:, 1], alpha=0.1)

forecast.predicted_mean.plot(ax=ax, style='r', label='Forecast')
forecast_ci = forecast.conf_int()
forecast_index = np.arange(len(predict_ci), len(predict_ci) + len(forecast_ci))
ax.fill_between(forecast_index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], alpha=0.1)

# Cleanup the image
ax.set_ylim((4, 8));
legend = ax.legend(loc='lower left');

References¶

Commandeur, Jacques J. F., and Siem Jan Koopman. 2007.
An Introduction to State Space Time Series Analysis.
Oxford ; New York: Oxford University Press.

Durbin, James, and Siem Jan Koopman. 2012.
Time Series Analysis by State Space Methods: Second Edition.
Oxford University Press.

最終更新日: 2025年01月28日