wwrechard / pydlm

A python library for Bayesian time series modeling
BSD 3-Clause "New" or "Revised" License
475 stars 98 forks source link

Univariate DLM Approach #14

Open tian-yu-ic opened 7 years ago

tian-yu-ic commented 7 years ago

Hi,

Amazing work here! I am reading Bayesian Forecasting and Dynamic Models (Harrison and West, 1997) and try to use your module to construct an univariate DLM model as follows. image It seems the components that I need to put into the model are just trends, dynamic and autoreg, but I am wondering how to mix them so that I can get such a model. Is there a way to do so? Because I can not find any similar example in the documentation.

Thanks!

tian-yu-ic commented 7 years ago

I think it may help to specify my model as follows. image This is actually what I want to construct.

wwrechard commented 7 years ago

Hi tian-yi-ic,

Per the case you provided, if I understand correctly, you want to build a model of 1-lag auto-regression + constant trend? You can definitely fit such model using PyDLM. Following is a simple fake example.

import numpy as np
from pydlm import dlm, autoReg, trend

# Generate fake random walk data using numpy.
n = 100
off_set = 0
y = [0] * n
y[0] = off_set + np.random.normal() * 2
for i in range(1, n):
  constant = off_set + np.random.normal() * 2
  y[i] = y[i - 1] + constant + np.random.normal()

# Build model with pydlm
myDLM = dlm(y) + trend(degree=0, discount=0.98) + autoReg(degree=1, data=y, discount=1.0)
myDLM.fit()
myDLM.plot()
tian-yu-ic commented 7 years ago

Thank you so much! That's exactly what I want.

tian-yu-ic commented 7 years ago

Sorry to bother you again for the same problem. I think image contributes to my model and should be included in the pydlm model. But why doesn't your example take it into account? Actually I am asking that the dlm model consists of observation and system(on the left side of the first graph), how does pydlm deal with the system part?

wwrechard commented 7 years ago

Phi/Gt is automatically generated in pydlm. In this special example, phi_1 and phi_2 will be multiplied to theta_1 and theta_2, so they are unidentifiable, i.e., (phi_1 = 1, theta_1 = 1) and (phi_1 = 0.5, theta_1 = 2) gives you the same result (similar for phi_2 and theta_2). Since theta_1 and theta_2 are the latent states and going to learnt by the DLM from the data, phi_1 and phi_2 will just be fixed to be 1 (because no matter what value you specified, it will be absorbed by theta_1 and theta_2 when fitting the model).

For more sophisticated cases such as polynomial trend and seasonality, Phi_1 and Phi_2 will take more complicated forms rather than an identity matrix like in this case. In general, the transition matrix for each type of component is fixed a-priori in DLM. PyDLM will automatically generate these fixed transition matrix and combine them to construct the system transition matrix Phi. You don't need to worry about it or supply it yourself. For better understanding of the transition matrix (Phi/Gt) and the canonical form, you can read the chapter 7, 8 and 9 of the book Bayesian Forecasting and Dynamic Models (Harrison and West, 1999)

In most cases you don't even need to worry about phi_1 and phi_2. If you are interested in knowing the actual values, you can reference to the code under pydlm/modeler/... and look into each component, or you can query lively from the constructed model

trend3 = trend(degree=3)
trend3.transition

If you do want to supply your own transition matrix (phi_1 and phi_2), you can directly modify when constructing the component.

phi_1 = 1
phi_2 = 1
trend0 = trend(degree=0, discount=0.98)
trend0.transition = np.matrix(phi_1)
ar1 = autoReg(degree=1, data=y, discount=1.0)
ar1.transition = np.matrix(phi_2)
myDLM = dlm(y) + trend0 + ar1
myDLM.fit()
myDLM.plot()

But I highly recommend not doing this without a good understanding of the transition matrix in DLM. The transition matrix must satisfy some conditions (e.g., spectral_norm = 1) to be eligible. You can try some random transition matrix and see their weird behavior in the model fitting.

For the latent states (theta's) you can get them easily after fitting the model

myDLM.getLatentState(filterType='forwardFilter')
myDLM.getLatentState(filterType='backwardSmoother')