Dynamic linear models — user manual¶
This package implements the Bayesian dynamic linear model (DLM, Harrison and West, 1999) for time series analysis. The DLM is built upon two layers. The first layer is the fitting algorithm. DLM adopts a modified Kalman filter with a unique discounting technique from Harrison and West (1999). Like the usual Kalman filter, it accepts a transition matrix, a measurement matrix, an observation, a latent state, an innovation and an error covariance matrix and return the updated state and error covariance. These quantities will all be supplied internally – users are free from any annoying calculations. Different from the usual Kalman filter, the modified Kalman filter does not require the tuning of the two parameters: the error covariance matrix and the observational variance, so the model fitting is extremely efficient (could be up to 1000 times faster than the EM algorithm), more details will be provided in the section of the discounting technique.
The second layer of DLM is itsmodeling feature. Nicely summarized in Harrison and West (1999), most common models can be expressed in one unified form – canonical form, which is closely related to the Jordan decomposition. Thanks to this keen observation, the DLM can easily incorporate most modeling components and turn them into the corresponding transition matrices and other quantities to be supplied to the Kalman filter. Examples are trend, seasonality, holidays, control variables and auto-regressive, which could appear simultaneously in one model. Due to this nice property, users of this package can construct models simply by “adding” some component into the model as:
myDLM = dlm(data) + trend(1)
The modeling process is simple.
The purpose of the modeling is to better understand the time series data and for to forecast into the future. So the key output from the model are the filtered time series, smoothed time series and one-step ahead prediction. We will cover this topic later in this section.
The advantage of pydlm:
flexibility in constructing complicated models
Extremely efficient model fitting with the discounting technique
user-specific adjustment on the adaptive property of the model
The disadvantage of pydlm:
only for Gaussian noise
Miscellaneous¶
PyDLM index starts at 0 instead of 1, i.e., for any prediction or modification that involves date argument, it corresponds to the actual index in the array. For instance, for a time series with length 10, the date of the last day is 9.
Modeling¶
As discussed in the beginning, the modeling process is very simple
with pydlm, most modeling functions are integrated in the class
dlm. Following is an example for constructing a dlm with
linear trend, 7-day seasonality and control variables:
from pydlm import dlm, trend, seasonality, dynamic, autoReg, longSeason
data = [0] * 100 + [3] * 100
SP500Index = [[2000] for i in range(100)] + [[2010] for i in range(100)]
page = [[i, i + 1, i + 2, i + 3] for i in range(200)]
myDLM = dlm(data)
myDLM = myDLM + trend(degree=1, discount=0.95, name='trend1')
myDLM = myDLM + seasonality(period=7, discount=0.99, name='week')
myDLM = myDLM + dynamic(features=SP500Index, discount=1, name='SP500')
myDLM = myDLM + dynamic(features=page, discount=1, name='page')
User can also use dlm.add() method to add new component:
myDLM.add(trend(degree=0, discount=0.99, name='trend2'))
The imput data must be an 1d array or a list, since the current
dlm only supports one dimensional time series. Supporting for
multivariate time series will be built upon this one dimensional class
and added in the future.
Now the variable myDLM contains the data and the modeling information. It will construct the corresponding transition, measurement, innovation, latent states and error covariance matrix once model fitting is called. Modify an existing model is also simple. User can brows the existing components of the model by:
myDLM.ls()
It will show all the existing components and their corresponding names. Name can be specified when the component is added to the dlm, for example:
myDLM = myDLM + seasonality(4, name = 'day4')
myDLM.ls()
We can also easily delete the unwanted component by using delete:
myDLM.delete('trend2')
myDLM.delete('day4')
After the building steps, we can specify some parameters for the model fitting, the most common one is the prior guess on the observational noise. The default value is 1.0. To change that to 10 you can do:
myDLM.noisePrior(10.0)
Such change usually has small impact on the model and is almost ignorable.
Model components¶
There are four model components provided with this package: trend, seasonality, dynamic and the auto-regression.
Trend¶
trend class is a model component for trend
behavior. The data might be const or increasing linearly or
quadraticly, which can all be captured by trend. The degree
argument specifics the shape of the trend. degree=0 indicates this
is a const, degree=1 indicates a line and degree=2 stands for a
quadratic curve and so on so forth. w sets the prior
covariance for the trend component (same for all the other
components). The discounting factor will be explained later in next
section:
linearTrend = trend(degree=1, discount=0.99, name='trend1', w=1e7)
Seasonality¶
The seasonality class models the periodic behavior of the
data. Compared to the sine or cosine periodic curves,
seasonality in this packages is more flexible, since it can
turn into any shapes, much broader than the triangular families:
weekPeriod = seasonality(period=7, discount=0.99, name='week', w=1e7)
In the package, we implements the seasonality component in a
form-free way (Harrison and West, 1999) to avoid the identifiability
issue. The states of one seasonality component are always summed up to
zero, so that it will not tangle with the trend component.
Dynamic¶
The dynamic class offers the modeling ability to add any additional
observed time series as a controlled variable to the current one. For
example, when studying the stock price, the ‘SP500’ index could be a
good indicator for the modeling stock. A dynamic component need the
user to supply the necessary information of the control variable over
time:
SP500 = dynamic(features=SP500Index, discount=0.99, name='SP500', w=1e7)
The input features for dynamic should be a list of
lists, since multi-dimension features are allowed. Following is one
simple example:
Features = [[2000], [2010], [2020], [2030]]
Features = [[1.0, 2.0], [1.0, 3.0], [3.0, 3.0]]
Auto-regression¶
The autoReg class constructs the auto-regressive component on
the model, i.e., the direct linear or non-linear dependency between
the current observation and the previous days. User needs to specify
the number of days of the dependency:
AR3 = autoReg(degree=3, discount=0.99, name='ar3', w=1e7)
There is once a data argument needed for constructing autoregression
features but is now deprecated. autoReg can now fetch the
data directly from the main dlm class and no need to provide
during instantiation.
In this example, the latent stats for Auto-regression are aligned in a way of [today - 3, today - 2, today - 1]. So when fetching the coefficients from the latent states, this will be the correct order to read the coefficients.
Long-seasonality¶
The longSeason class is a complement class for
seasonality. It allows constructing seasonality component that
does not change every step. For example, the time unit is day, but
user wants to add a monthly seasonality, then longSeason is
the correct choice:
monthly = longSeason(period=12, stay=30, data=data, name='monthly', w=1e7)
These five classes of model components offer abundant modeling possiblities of the Bayesian dynamic linear model. Users can construct very complicated models using these components, such as hourly, weekly or monthly periodicy and holiday indicator and many other features.
Model fitting¶
Entailed before, the fitting of the dlm is fulfilled by a modified
Kalman filter. Once the user finished constructing the model by adding
different components. the dlm will compute all the necessary
quantities internally for using Kalman filter. So users can simply
call dlm.fitForwardFilter(), dlm.fitBackwardSmoother() or
even simply dlm.fit() to fit both forward filter and backward
smoother:
myDLM.fitForwardFilter()
myDLM.fitBackwardSmoother()
myDLM.fit()
The dlm.fitForwardFilter() is implemented in an online
manner. It keeps an internal count on the filtered dates and once new
data comes in, it only filters the new data without touching the
existing results. In addition, this function also allows a rolling
window fitting on the data, i.e., there will be a moving window and
for each date, the kalman filter will only use the data within the
window to filter the observation. This is equivalent to that the model
only remembers a fixed length of dates:
myDLM.fitForwardFilter(useRollingWindow=True, windowLength=30)
myDLM.fitBackwardSmoother()
For dlm.backwardSmoother(), it has to use the whole time series
to smooth the latent states once new data comes in. The smoothing
provides a good retrospective analysis on our past decision of the
data. For example, we might initially believe the time series is
stable, while that could be a random behavior within a volatile time
series, and the user learn this from the smoother.
Once the model fitting is completed, users can fetch the filtered or
smoothed results from dlm:
myDLM.getMean(filterType='forwardFilter')
myDLM.getMean(filterType='backwardSmoother')
myDLM.getMean(filterType='predict')
myDLM.getVar(filterType='forwardFilter')
myDLM.getVar(filterType='backwardSmoother')
myDLM.getVar(filterType='predict')
The dlm recomputes a wide variety of model quantities that
can be extracted by the user. For example, user can get the filtered
states and covariance by typing:
myDLM.getLatentState(filterType='forwardFilter')
myDLM.getLatentState(filterType='backwardSmoother')
myDLM.getLatentCov(filterType='forwardFilter')
myDLM.getLatentCov(filterType='backwardSmoother')
This can be specified into individual component. For example, assume the
model contains a trend component with a name of trend1, we
can extract the corresponding latent state only for trend1 as:
myDLM.getLatentState(filterType='forwardFilter', name='trend1')
myDLM.getLatentState(filterType='backwardSmoother', name='trend1')
myDLM.getLatentCov(filterType='forwardFilter', name='trend1')
myDLM.getLatentCov(filterType='backwardSmoother', name='trend1')
as well as the mean of trend1 (evaluation * latent states):
myDLM.getMean(filterType='forwardFilter', name='trend1')
myDLM.getVar(filterType='forwardFilter', name='trend1')
One can also get the confidence interval on the filtered time series:
myDLM.geInterval(filterType='forwardFilter', p = 0.99)
There are also corresponding methods for smoothed and predicted
results. For more detail, please refer to the dlm class
documentation.
Model prediction¶
dlm provides three predict functions: dlm.predict() and
dlm.continuePredict() and dlm.predictN(). The last one is
wrapper of the former two and is recommended to use. (The former two
will be gradually deprecated).
The dlm.predict() is a one-day ahead
prediction function based on a user given date and feature set:
# predict next date after the time series
featureDict = {'SP500':[2090], 'page':[1, 2, 3, 4]}
(predictMean, predictVar) = myDLM.predict(date=myDLM.n - 1, featureDict=featureDict)
The function returns a tuple of predicted mean and predicted variance.
The featureDict argument is a dictionary contains the feature
information for dynamic component. Suppose the model contains
a one-dimensional dynamic component named SP500 and another
four-dimensional dynamic component page, then the featureDict takes the
following Form:
featureDict = {'SP500':[2090], 'page':[1, 2, 3, 4]}
If the featureDict is not supplied but the date is not the last day, then the algorithm will automatically fetch from the old data about the feature value of all the dynamic component:
# predict a day in the middle
(predictMean, predictVar) = myDLM.predict(date=myDLM.n - 10)
The algorithm will use the feature on the date of myDLM.n - 9 in featureDict. If date is the last day but the featureDict is not provided, then an error will be raised.
If the user is interested beyond one-day ahead prediction, they can
use dlm.continuePredict() for multiple-day ahead prediction,
after using dlm.predict():
feature1 = {'SP500':[2090], 'page':[10, 20, 30, 40]}
feature2 = {'SP500':[2010], 'page':[11, 21, 31, 41]}
feature3 = {'SP500':[1990], 'page':[12, 22, 32, 42]}
# one-day ahead prediction after the last day
(predictMean, predictVar) = myDLM.predict(date=myDLM.n - 1, featureDict=feature1)
# we continue to two-day ahead prediction after the last day
(predictMean, predictVar) = myDLM.continuePredict(featureDict=feature2)
# we continue to three-day ahead prediction after the last day
(predictMean, predictVar) = myDLM.continuePredict(featureDict=feature3)
dlm.continuePredict() can only be used after dlm.predict()
for multiple-day prediction. The featureDict can also be ignored if
the prediction is requested on dates before the last day and the
features on the predict day can be found from the old data.
dlm.predictN() (recommended) predicts over multiple days and is
a wrapper of the two functions above. Using the same example, the
results can be obtained by just called:
features = {'SP500':[[2090], [2010], [1990]], 'page':[[10, 20, 30,
40], [11, 21, 31, 41], [12, 22, 32, 42]]}
(predictMean, predictVar) = myDLM.predictN(N=3, date=myDLM.n-1,
featureDict=features)
The predictMean and predictVar will be two lists of three elements containing the predicted mean and variance.
Model amending¶
The user can still add, delete, modify data even when the model has been constructed.
Adding new data¶
For adding more data, user can opt to
dlm.append():
newData = [0, 1, 2]
myDLM.append(newData, component = 'main')
If the model contains dynamic component, the corresponding
features need to be updated as well:
newSP500 = [[2000], [2100], [2200]]
myDLM.append(data = newSP500, component = 'SP500')
Then the user can rerun the forward filter:
myDLM.fitForwardFilter()
The package will continue running the forward filter on the three new data ponts.
Deleting existing data¶
To delete any existing data, user can simply use the dlm.popout()
function from dlm on a specific date, for example:
myDLM.popout(0)
Different from dlm.append(), dlm.popout() will be executed
automatically for all components, so the user does not need to conduct
the deletion mannually for each component. After the deletion, the
forward filter needs to be rerun following the deleted date:
myDLM.fitForwardFilter()
Again, the package will automatically recognize the date and fit only the necessary period of time.
Ignoring a date¶
Ignoring is very similar to deleting. The only difference is the time
counts. Because deleting will delete the data entirely, the time
counts will therefore reduce by 1. By contrast, ignoring will treat
the specific date as missing data, so the time count will not
change. This difference is important when preriodicy is
concerned. Changing of time counts will have high impacts on
seasonality components.
dlm.ignore() simply set the data of a specific date to be None:
myDLM.ignore(2)
modify data¶
The dlm also provides user the ability to modify the data on a
specific date and a specific component. This function enables possible
future extension to interactive anomaly detection and data debugging:
myDLM.alter(date = 2, data = 0, component = 'main')
Model plotting¶
This package offers rich ploting options for illustrating the
results. User can simply call dlm.plot() to directly plot the
results once the models are fitted:
myDLM.plot()
# plot the mean of a given component
myDLM.plot(name=the_component_name)
# plot the latent state of a given component
myDLM.plotCoef(name=the_component_name)
User can choose which results to plot via dlm.turnOn() and
dlm.turnOff():
myDLM.turnOn('filtered plot')
myDLM.turnOff('predict plot')
myDLM.turnOff('smoothed plot')
User can also choose whether to plot the confidence interval and whether plot the results in one figure or separate figures. The default is to plot the confidence interval and in separate plots. To change that:
myDLM.turnOff('confidence')
myDLM.turnOff('multiple plots')
The quantile of the confidence interval can be set via
dlm.setConfidence():
myDLM.setConfidence(p = 0.95)
Currently there are two types of confidence interval realization. The default is ‘ribbon’ and the alternative is ‘line’. Users can change the confidence interval plots by:
myDLM.setIntervalType('ribbon')
myDLM.setIntervalType('line')
The default colors for the plots are:
original data: ‘black’
filtered results: ‘blue’
one-step ahead prediction: ‘green’
smoothed results: ‘red’
User can change the color setting via dlm.setColor(). The color
space is the same as the matplotlib:
myDLM.setColor('filtered plot', 'yellow')
myDLM.setColor('data', 'blue')
If user decide to go back to the original setting, they can use
dlm.resetPlotOptions() to reset all the plot option:
myDLM.resetPlotOptions()
Model Tuning and evaluation¶
The discounting factor of DLM determines how fast the model adapts to
the new data. It could cause troublesome when users actually want the
model itself to figure that out. The modelTuner provides
users automatic tool for tuning the discounting factors:
from pydlm import modelTuner
myTuner = modelTuner(method='gradient_descent', loss='mse')
tunedDLM = myTuner.tune(myDLM, maxit=100)
The tuned discounting factor will be set inside the tunedDLM. Users can examine the tuned value from myTuner:
tuned_discounts = myTuner.getDiscounts()
After tuning, myDLM will remain unchanged and tunedDLM will contain the tuned discounting factos and will be in uninitialized status. Users need to run:
tunedDLM.fit()
before any other analysis on tunedDLM. The dlm also provides
a simpler way to tune the discounting factor if the user would like
myDLM to be altered directly:
myDLM.tune()
The tuner makes use of the MSE (one-day ahead prediction loss) and the gradient descent algorithm to tune the discounting factor to achieve the minimum loss. The discounting factors are assumed to be different across components but the same within a component. For now, only the MSE loss and the gradient descent algorithm is supported.
To ease the evaluation of the performance of the model fitting, the model also provides the residual time series and the one-day a head prediction error via:
mse = myDLM.getMSE()
residual = myDLM.getResidual(filterType='backwardSmoother')
Users can use these values to evaluate and choose the optimal model.
Advanced Settings¶
This part of settings closely relate to the algorithm behavior and
offers some advanced features, some of which are still under
development. Currently implemented is the dlm.stableMode()
function and the dlm.evolveMode(). The dlm.stableMode() is
turned on by default, and you can turn it off by:
myDLM.stableModel(False)
This mode helps increasing the numerical stability of the dlm
when small discounting factor is used. Details about discounting
factor will be covered in next section. The dlm.evolveMode() is
used to control how different components evolve over time. See
Harrison and West (1999, Page 202). They could
evolve independently, which is equivalent to assume the innovation
matrix is a block-diagonal matrix. The default assumption is
‘independent’ and to change to ‘dependent’, we can simply type:
myDLM.evolveMode('dependent')
The difference between ‘independent’ and ‘dependent’ is best explained when there are multiple components with different discounting factor and one of them is One. In the ‘dependent’ mode, the smoothed latent states of the component with discount factor 1 will be a value fluctuating around a constant, while in the ‘independent’ mode, it will be an exact constant. User can choose which to use depending on their own use case.
In the future, following functionalities are planned to be added: feature selection among dynamic components, factor models for high dimensional latent states.