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.3305405227418
 1132.7806610433695
 1068.0283899362175
 1113.9581805683022
 1127.584244381744
 1136.7260100000785
 1047.7788084559295
 1097.1469834323289
 1170.5115745331705
 1162.3410382738443
    ⋮
  919.1502104663157
  914.3128939030667
  982.4575292298758
  963.6794895051129
  905.6644817994512
  909.218607146413
  858.2558407769868
  819.8093902947194
  798.5389325062581

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.1898201993902
 1107.569196194333
 1102.871214871117
 1111.7182506919808
 1111.047760118626
 1105.6367735894357
 1094.9612954747747
 1111.5907248401259
 1116.7534942244624
 1097.3920735131546
    ⋮
  914.7071002120507
  913.0926705484638
  912.649295374962
  887.2840721307691
  859.5253578713462
  842.7604418420963
  818.6124852722882
  804.2078440573923
  798.5389325062581

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.5389325062581
 798.5389325062581
 798.5389325062581
 798.5389325062581
 798.5389325062581
 798.5389325062581
 798.5389325062581
 798.5389325062581
 798.5389325062581
 798.5389325062581

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.8126525851467
 1098.7561144298797
 1096.3829253092656
 1099.0551652923414
 1097.5307339500566
 1093.643299605073
 1087.2458136427542
 1091.2221278239701
 1089.948940139774
 1078.08235930078
    ⋮
  896.679952627532
  894.2133507269256
  891.4900326261046
  878.2316157351054
  863.695854023848
  853.6121320762915
  841.0550061972679
  833.1526388503645
  829.7574165555225

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.1291     -3.1079       0.0019
s_ma_L12        -0.5561         0.1462     -3.8029       0.0001
sigma2_η         0.0014         0.0024      0.5646       0.5724

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

forec = forecast(model, 24)

StateSpaceModels.Forecast{Float64}([[6.1100783269760806], [6.053643535089317], [6.171326207071635], [6.199284566888648], [6.232663436588682], [6.368794427137307], [6.50742299366151], [6.503004703603158], [6.324632038853844], [6.209032320028253]  …  [6.26744855378555], [6.295406913602562], [6.328785783302597], [6.464916773851222], [6.603545340375425], [6.5991270503170725], [6.420754385567759], [6.305154666742168], [6.159527964824452], [6.263985795085841]], [[0.0013606492643246676;;], [0.0018485279562723811;;], [0.0023364066482185275;;], [0.002824285340164674;;], [0.0033121640321108203;;], [0.0038000427240569667;;], [0.004287921530288826;;], [0.004775800322234974;;], [0.0052636791141811206;;], [0.005751557906128267;;]  …  [0.010255711840981315;;], [0.0112756523907784;;], [0.012295592969832864;;], [0.013315533581722543;;], [0.01433547439154342;;], [0.015355415366095434;;], [0.01637535614434456;;], [0.017395296922618417;;], [0.01842137626486757;;], [0.019434891531997466;;]])

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.006970992671995191
 -0.018017085598208898
 -0.05779311131016858
 -0.08804898191193647
 -0.08672338383008846
 -0.09885574748987086
 -0.10336477068885125
 -0.0911911269098084
 -0.06401598510954876
 -0.05203595297846356
  ⋮
 -0.044136693650863305
 -0.016113957869339043
 -0.00825586270775562
  0.0021714599595704376
 -0.0017053976021500876
 -0.001980639561218401
 -0.020603423812991215
 -0.03567384008771211
 -0.035703311976724876

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.034852443288206914 0.06796954095741636 … 0.020043103647071625 0.043596386754162; 0.08528143585934789 0.0605080674020293 … 0.057683912071698096 0.0049623533719485025; … ; 0.042618122885714094 0.15619450040736726 … 0.009869444793378257 0.057266585129702285; 0.12340374360452433 0.1543141909082557 … 0.039794928573065214 0.024643478514379424], [0.03134192170231125, 0.040697685364027854, 0.047771990949174115, 0.050466803719844855, 0.05194281017497467, 0.054719743531486426, 0.05634143632507803, 0.05677756399576602, 0.05769456750211524, 0.05794235840320962  …  0.07949400164891586, 0.08412011378980115, 0.08789479633685823, 0.09059329873866144, 0.09059711681065326, 0.08991636041781906, 0.0907243129559548, 0.08925213241730257, 0.08571336825122462, 0.08986521990849447], [0.021760345428662963 0.05072895075113007 … 0.011904030294759376 0.02922773673975335; 0.06648950106458322 0.042751562773980935 … 0.03786198842044632 0.008900692647500383; … ; 0.02788496727404974 0.11289420376397297 … 0.027534829085642264 0.03868424087803118; 0.08395030769302766 0.10971736327969013 … 0.033142958993894585 0.02988381379588632], [0.02213533160563259, 0.029153349827080094, 0.03499671782466187, 0.037863548338893084, 0.03949557130610761, 0.041289799690621236, 0.04229558760058205, 0.043692281006359884, 0.04390206485589362, 0.04469632818152791  …  0.05732320471459466, 0.06122944118526715, 0.06335240076717749, 0.06516450879603486, 0.06608409535305212, 0.06503160316624837, 0.06533840949145603, 0.06419038993128656, 0.06215454001593312, 0.06290885224344045])