Examples

In this section we show examples of applications and use cases of the package.

Nile river annual flow

Here we will follow an example from Durbin & Koopman's book. We will use the LocalLevel model applied to the annual flow of the Nile river at the city of Aswan between 1871 and 1970.

First, we load the data:

using CSV, DataFrames
nile = CSV.File(StateSpaceModels.NILE) |> DataFrame

Next, we fit a LocalLevel model:

model = LocalLevel(nile.flow)
fit!(model)

LocalLevel

We can analyze the filtered estimates for the level of the annual flow:

filter_output = kalman_filter(model)
get_filtered_state(filter_output)

100×1 Matrix{Float64}:
 1103.3305405019873
 1132.7806610192085
 1068.028390006868
 1113.958180514016
 1127.5842443027625
 1136.7260099107473
 1047.7788087441866
 1097.1469834335082
 1170.5115742111407
 1162.3410380748621
    ⋮
  919.1502103715065
  914.3128938560766
  982.4575288777272
  963.6794893343584
  905.6644819446645
  909.2186072363559
  858.2558410805408
  819.8093906966042
  798.5389329001946

We can do the same for the smoothed estimates for the level of the annual flow:

smoother_output = kalman_smoother(model)
get_smoothed_state(smoother_output)

100×1 Matrix{Float64}:
 1107.1898201672168
 1107.569196156864
 1102.871214880668
 1111.7182506002334
 1111.0477600215625
 1105.636773529889
 1094.961295504984
 1111.5907246621625
 1116.753493947436
 1097.3920733804534
    ⋮
  914.7070999995475
  913.0926703034565
  912.649295060822
  887.2840719916562
  859.5253579201594
  842.7604419624364
  818.612485557153
  804.2078444270194
  798.5389329001946

StateSpaceModels.jl can also be used to obtain forecasts. Here we forecast 10 steps ahead:

steps_ahead = 10
forec = forecast(model, steps_ahead)
expected_value = forecast_expected_value(forec)

10×1 Matrix{Float64}:
 798.5389329001946
 798.5389329001946
 798.5389329001946
 798.5389329001946
 798.5389329001946
 798.5389329001946
 798.5389329001946
 798.5389329001946
 798.5389329001946
 798.5389329001946

We can also simulate multiple scenarios for the forecasting horizon based on the estimated predictive distributions.

scenarios = simulate_scenarios(model, 10, 100)
size(scenarios)

(10, 1, 100)

The package also handles missing values automatically. To that end, the package considers that any NaN entries in the observations are missing values.

nile.flow[[collect(21:40); collect(61:80)]] .= NaN

40-element view(::Vector{Float64}, [21, 22, 23, 24, 25, 26, 27, 28, 29, 30  …  71, 72, 73, 74, 75, 76, 77, 78, 79, 80]) with eltype Float64:
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
   ⋮
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN
 NaN

Even though the series has several missing values, the same analysis is possible:

model = LocalLevel(nile.flow)
fit!(model)

LocalLevel

And the exact same code can be used for filtering and smoothing:

filter_output = kalman_filter(model)
smoother_output = kalman_smoother(model)
get_smoothed_state(smoother_output)

100×1 Matrix{Float64}:
 1098.8126525445964
 1098.7561143909682
 1096.3829252751582
 1099.0551652516476
 1097.5307339096387
 1093.6432995681325
 1087.2458136129192
 1091.2221277793437
 1089.9489400891773
 1078.0823592635072
    ⋮
  896.679952593925
  894.2133506941585
  891.4900325934517
  878.2316157224227
  863.6958540332128
  853.6121320991451
  841.0550062394265
  833.1526389041042
  829.7574166138527

Airline passengers

We often write the model SARIMA model as an ARIMA $(p,d,q) \times (P,D,Q,s)$, where the lowercase letters indicate the specification for the non-seasonal component, and the uppercase letters indicate the specification for the seasonal component; $s$ is the periodicity of the seasons (e.g. it is often 4 for quarterly data or 12 for monthly data). The data process can be written generically as

\[\begin{equation} \phi_p (L) \tilde \phi_P (L^s) \Delta^d \Delta_s^D y_t = A(t) + \theta_q (L) \tilde \theta_Q (L^s) \epsilon_t \end{equation}\]

where

• $\phi_p (L)$ is the non-seasonal autoregressive lag polynomial,
• $\tilde \phi_P (L^s)$ is the seasonal autoregressive lag polynomial,
• $\Delta^d \Delta_s^D y_t$ is the time series, differenced $d$ times, and seasonally differenced $D$ times.,
• $A(t)$ is the trend polynomial (including the intercept),
• $\theta_q (L)$ is the non-seasonal moving average lag polynomial,
• $\tilde \theta_Q (L^s)$ is the seasonal moving average lag polynomial

sometimes we rewrite this as:

\[\begin{equation} \phi_p (L) \tilde \phi_P (L^s) y_t^* = A(t) + \theta_q (L) \tilde \theta_Q (L^s) \epsilon_t \end{equation}\]

where $y_t^* = \Delta^d \Delta_s^D y_t$. This emphasizes that just as in the simple case, after we take differences (here both non-seasonal and seasonal) to make the data stationary, the resulting model is just an ARMA model.

As an example, consider the airline model ARIMA $(0,1,1) \times (0,1,1,12)$. The data process can be written in the form above as:,

\[\begin{equation} \Delta \Delta_{12} y_t = (1 - \theta_1 L) (1 - \tilde \theta_1 L^{12}) \epsilon_t \end{equation}\]

Here, we have:

• $\phi_p (L) = 1$, (i.e. there is no auto regressive effect)
• $\tilde \phi_P (L^s) = 1$, (i.e. there is no seasonal auto regressive effect)
• $d = 1, D = 1, s=12$ indicating that $y_t^*$ is derived from $y_t$ by taking first-differences and then taking 12-th differences.,
• $A(t) = 0$ no trend,
• $\theta_q (L) = (1 - \theta_1 L)$,
• $\tilde \theta_Q (L^s) = (1 - \tilde \theta_1 L^{12})$,

It may still be confusing to see the two lag polynomials in front of the error variable, but notice that we can multiply the lag polynomials together to get the following model:

\[\begin{equation} y_t^* = (1 - \theta_1 L - \tilde \theta_1 L^{12} + \theta_1 \tilde \theta_1 L^{13}) \epsilon_t \end{equation}\]

For the airline model ARIMA $(0,1,1) \times (0,1,1,12)$ with an intercept, the command is:

using CSV, DataFrames
air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
log_air_passengers = log.(air_passengers.passengers)
model = SARIMA(log_air_passengers; order = (0, 1, 1), seasonal_order = (0, 1, 1, 12))
fit!(model)
print_results(model)

                             Results
===============================================================
Model:                        SARIMA(0, 1, 1)x(0, 1, 1, 12) with zero mean
Number of observations:       144
Number of unknown parameters: 3
Log-likelihood:               244.6323
AIC:                          -483.2645
AICc:                         -483.0931
BIC:                          -474.3551
---------------------------------------------------------------
Parameter      Estimate      Std.Error      z stat      p-value
ma_L1           -0.4011         0.1178     -3.4057       0.0007
s_ma_L12        -0.5561         0.1306     -4.2593       0.0000
sigma2_η         0.0014         0.0024      0.5655       0.5717

To make a forecast of 24 steps ahead of the model the command is:

forec = forecast(model, 24)

StateSpaceModels.Forecast{Float64}([[6.110078764455084], [6.053641814899042], [6.171319394066974], [6.199285353951767], [6.232664971514262], [6.368794028735144], [6.5074254773101785], [6.503005891849219], [6.324630808423989], [6.209033497289308]  …  [6.267440442347349], [6.295406402232142], [6.328786019794636], [6.464915077015519], [6.603546525590553], [6.599126940129594], [6.4207518567043635], [6.305154545569682], [6.159525754508261], [6.2639826897325905]], [[0.0013606604321276265;;], [0.0018485288142440796;;], [0.002336397196358967;;], [0.0028242655784738537;;], [0.0033121339605887406;;], [0.0038000023427036274;;], [0.0042878708391042245;;], [0.004775739321219113;;], [0.005263607803334;;], [0.005751476285449886;;]  …  [0.01025591850196266;;], [0.011275928751177399;;], [0.012295939029649608;;], [0.013315949340957433;;], [0.01433595985019773;;], [0.015355970524169855;;], [0.016375981001836615;;], [0.017395991479528108;;], [0.018422136741535772;;], [0.019435726686224074;;]])

The text from this example is based on Python`s statsmodels library. The estimates for this example match up to the 3th decimal place the results of the paper State Space Methods in Ox/SsfPack from the journal of statistical software.

Finland road traffic fatalities

In this example, we will follow what is illustrated on Commandeur, Jacques J.F. & Koopman, Siem Jan, 2007. "An Introduction to State Space Time Series Analysis," OUP Catalogue, Oxford University Press (Chapter 3). We will study the LocalLinearTrend model with a series of the log of road traffic fatalities in Finalnd and analyse its slope to tell if the trend of fatalities was increasing or decrasing during different periods of time.

using StateSpaceModels, CSV, DataFrames
df = CSV.File(StateSpaceModels.VEHICLE_FATALITIES) |> DataFrame
log_ff = log.(df.ff)

34-element Vector{Float64}:
 6.961296045910167
 7.04141166379481
 7.052721049232323
 6.990256500493881
 6.762729506931879
 6.813444599510896
 6.6895992691789665
 6.5638555265321274
 6.413458957167357
 6.476972362889683
 ⋮
 6.089044875446846
 6.0014148779611505
 6.082218910376446
 5.991464547107982
 6.066108090103747
 5.981414211254481
 6.07073772800249
 6.028278520230698
 5.937536205082426

We fit a LocalLinearTrend

model = LocalLinearTrend(log_ff)
fit!(model)

LocalLinearTrend

By extracting the smoothed slope we conclude that according to our model the trend of fatalities in Finland was increasing in the years 1970, 1982, 1984 through to 1988, and in 1998

smoother_output = kalman_smoother(model)
get_smoothed_state(smoother_output)[:, 2]

34-element Vector{Float64}:
  0.006970445554088656
 -0.018017366528511047
 -0.05779302153277158
 -0.08804871096699536
 -0.08672341922134437
 -0.09885570021686363
 -0.10336462347585701
 -0.09119102208103061
 -0.06401607497924203
 -0.0520358994854154
  ⋮
 -0.04413673123532372
 -0.016114139111477396
 -0.008255914333836822
  0.0021713539702367148
 -0.0017054233098589694
 -0.0019807509373179436
 -0.02060331834406779
 -0.03567353864944685
 -0.035703010022977676

Cross validation of the forecasts of a model

Often times users would like to compare the forecasting skill of different models. The function cross_validation makes it easy to make a rolling window scheme of estimations and forecasts that allow users to track each model forecasting skill per lead time.

using CSV, DataFrames
air_passengers = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
log_air_passengers = log.(air_passengers.passengers)
model = BasicStructural(log_air_passengers, 12)
b = cross_validation(model, 24, 50)

StateSpaceModels.CrossValidation{Float64}([0.03484377234367386 0.06796426432653924 … 0.020043255441722252 0.043594053520111586; 0.08527276552221164 0.060501029258772654 … 0.057685327172803724 0.004961051963599417; … ; 0.04263047044630941 0.1562286461706428 … 0.009865862864648633 0.0572689454667481; 0.1234160920120555 0.15435007522521005 … 0.039793629051112944 0.02464599616377683], [0.030957445584884093, 0.04103408249217823, 0.048014036368352046, 0.05096733929791111, 0.05253689865920562, 0.0553087305123248, 0.05688724031259575, 0.05664759035012117, 0.057419986990728374, 0.05711830164114173  …  0.07853446480441988, 0.08322704110071356, 0.08687089226199164, 0.08949980503142439, 0.08981712193757098, 0.09101741399409163, 0.08931297444727791, 0.08782661668242861, 0.08457644050837021, 0.08851573904138911], [0.02183822098139513 0.05122902437959376 … 0.012210759255072326 0.028760682144761976; 0.06659123794247085 0.04269122657100135 … 0.03862349140373144 0.009222024857850782; … ; 0.02785500819915169 0.11343166155540747 … 0.02747808927397871 0.037397251924312216; 0.08201647074966673 0.11078520585934679 … 0.033467804051313006 0.029607723404283436], [0.02186968196935341, 0.029453851015919255, 0.035149154333483354, 0.03828408731330641, 0.039900108853270964, 0.041477168197050655, 0.04246587669760547, 0.04357666069112843, 0.043857966421140177, 0.04438662834999094  …  0.056987868864365916, 0.06093353203265952, 0.06294249880870352, 0.06456758498755119, 0.0656734823119057, 0.06533733250776659, 0.06465702433003868, 0.06354856552014894, 0.06153984585694406, 0.062001720845728366])