StateSpaceModels.jl

Build Status	Coverage	Documentation

StateSpaceModels.jl is a package for modeling, forecasting, and simulating time series in a state-space framework. Implementations were made based on the book "Time Series Analysis by State Space Methods" (2012) by James Durbin and Siem Jan Koopman. The notation of the variables in the code also follows the book.

Quickstart

import Pkg

Pkg.add("StateSpaceModels")

using StateSpaceModels

y = randn(100)

model = LocalLevel(y)

fit!(model)

print_results(model)

forecast(model, 10)

kf = kalman_filter(model)

v = get_innovations(kf)

ks = kalman_smoother(model)

alpha = get_smoothed_state(ks)

Features

Current features include:

Kalman filter and smoother
Maximum likelihood estimation
Forecasting and Monte Carlo simulation
User-defined models (user specifies the state-space system)
Several predefined models, including:
- Exponential Smoothing (ETS, all the linear ones)
- Unobserved components (local level, basic structural, ...)
- SARIMA
- Linear regression
- Naive models
Completion of missing values
Diagnostics for the residuals of fitted models
Visualization recipes

Quick Examples

Fitting and forecasting

Quick example of different models fit and forecast for the air passengers time-series

using CSV
using DataFrames
using Plots
using StateSpaceModels

airp = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
log_air_passengers = log.(airp.passengers)
steps_ahead = 30

# SARIMA
model_sarima = SARIMA(log_air_passengers; order = (0, 1, 1), seasonal_order = (0, 1, 1, 12))
fit!(model_sarima)
forec_sarima = forecast(model_sarima, steps_ahead)

# Unobserved Components
model_uc = UnobservedComponents(log_air_passengers; trend = "local linear trend", seasonal = "stochastic 12")
fit!(model_uc)
forec_uc = forecast(model_uc, steps_ahead)

# Exponential Smoothing
model_ets = ExponentialSmoothing(log_air_passengers; trend = true, seasonal = 12)
fit!(model_ets)
forec_ets = forecast(model_ets, steps_ahead)

# Naive model
model_naive = SeasonalNaive(log_air_passengers, 12)
fit!(model_naive)
forec_naive = forecast(model_naive, steps_ahead)

plt_sarima = plot(model_sarima, forec_sarima; title = "SARIMA", label = "");
plt_uc = plot(model_uc, forec_uc; title = "Unobserved components", label = "");
plt_ets = plot(model_ets, forec_ets; title = "Exponential smoothing", label = "");
plt_naive = plot(model_ets, forec_naive; title = "Seasonal Naive", label = "");

plot(plt_sarima, plt_uc, plt_ets, plt_naive; layout = (2, 2), size = (500, 500))

Automatic forecasting

Quick examples on automatic forecasting. When performing automatic forecasting users should provide the seasonal period if there is one.

model = auto_ets(log_air_passengers; seasonal = 12)
model = auto_arima(log_air_passengers; seasonal = 12)

Contributing

PRs such as adding new models and fixing bugs are very welcome!
For nontrivial changes, you'll probably want to first discuss the changes via issue.

Citing StateSpaceModels.jl

If you use StateSpaceModels.jl in your work, we kindly ask you to cite the following paper:

@article{SaavedraBodinSouto2019,
  title={StateSpaceModels.jl: a Julia Package for Time-Series Analysis in a State-Space Framework},
  author={Raphael Saavedra and Guilherme Bodin and Mario Souto},
  journal={arXiv preprint arXiv:1908.01757},
  year={2019}
}

Function to update just the Kalman filter

Maybe it is good to have a function to just update the Kalman filter. My sketch:

function update_kalman_filter(model::StateSpaceModel{Typ}, filter::FilterOutput, model_test::StateSpaceModel{Typ}; tol::Typ = Typ(1e-5)) where Typ
    time_invariant = model.mode == "time-invariant"

    # Predictive state and its covariance
    a = cat(filter.a, Matrix{Typ}(undef, model_test.dim.n, model.dim.m), dims=1)
    P = cat(filter.P, Array{Typ, 3}(undef, model.dim.m, model.dim.m, model_test.dim.n+1), dims=3)
    att = cat(filter.att, Matrix{Typ}(undef, model_test.dim.n, model.dim.m), dims=1)
    Ptt = cat(filter.Ptt, Array{Typ, 3}(undef, model.dim.m, model.dim.m, model_test.dim.n), dims=3)

    # Innovation and its covariance
    v = cat(filter.v, Matrix{Typ}(undef, model_test.dim.n, model.dim.p), dims=1)
    F = cat(filter.F, Array{Typ, 3}(undef, model.dim.p, model.dim.p, model_test.dim.n) , dims=3)

    # Kalman gain
    K = Array{Typ, 3}(undef, model.dim.m, model.dim.p, model.dim.n + model_test.dim.n)
    
    # Auxiliary structures
    ZP = Array{Typ, 2}(undef, model.dim.p, model.dim.m)
    P_Ztransp_invF = Array{Typ, 2}(undef, model.dim.m, model.dim.p)
    RQR = Array{Typ, 2}(undef, model.dim.m, model.dim.m)

    # Steady state initialization
    steadystate = filter.steadystate
    tsteady     = filter.tsteady

    # Initial state: big Kappa initialization
    StateSpaceModels.fill_a1!(a)
    StateSpaceModels.fill_P1!(P; bigkappa = Typ(1e6))
    y = vcat(model.y, model_test.y)
    d = vcat(model.d, model_test.d)
    Z = cat(model.Z, model_test.Z,dims=3)
    c = vcat(model.c, model_test.c)
    
    StateSpaceModels.update_ZP!(ZP, Z, P, filter.tsteady) # ZP = Z[:, :, t] * P[:, :, t]
    StateSpaceModels.update_F!(F, ZP, Z, model.H, filter.tsteady) # F_t = Z_t * P_t * Z_t + H
    StateSpaceModels.update_P_Ztransp_Finv!(P_Ztransp_invF, ZP, F, filter.tsteady) # P_Ztransp_invF   = ZP' * invF(F, t)
    StateSpaceModels.update_K!(K, P_Ztransp_invF, model.T, filter.tsteady) # K_t = T * P_t * Z_t * F^-1_t
    
    for t_ = 1:model_test.dim.n
        t = model.dim.n+t_-1
        if t in model_test.missing_observations
            steadystate  = false
            v[t, :]      = fill(NaN, model.dim.p)
            F[:, :, t]   = fill(NaN, (model.dim.p, model.dim.p))
            StateSpaceModels.update_att!(att, a, t) # attt = a_t
            StateSpaceModels.update_Ptt!(Ptt, P, t) # Pttt = P_t
            StateSpaceModels.update_a!(a, att, model.T, c, t) # a_t+1 = T * attt + c_t
            StateSpaceModels.update_P!(P, model.T, Ptt, RQR, t) # P_t+1 = T * Pttt * T' + RQR'
        elseif steadystate
            StateSpaceModels.update_v!(v, y, Z, d, a, t) # v_t = y_t - Z_t * a_t - d_t
            StateSpaceModels.repeat_matrix_t_plus_1!(F, t-1) # F[:, :, t]   = F[:, :, t-1]
            StateSpaceModels.repeat_matrix_t_plus_1!(K, t-1) # K[:, :, t]   = K[:, :, t-1]
            StateSpaceModels.update_att!(att, a, P_Ztransp_invF, v, t) # attt = a_t + P_t * Z_t * F^-1_t * v_t
            StateSpaceModels.repeat_matrix_t_plus_1!(Ptt, t-1) # Pttt = Pttt-1
            StateSpaceModels.update_a!(a, att, model.T, c, t) # a_t+1 = T * attt + c_t
            StateSpaceModels.repeat_matrix_t_plus_1!(P, t) # P__+1 = P_t
        else
            StateSpaceModels.update_v!(v, y, Z, d, a, t) # v_t = y_t - Z_t * a_t - d_t
            StateSpaceModels.update_ZP!(ZP, Z, P, t) # ZP = Z[:, :, t] * P[:, :, t]
            StateSpaceModels.update_F!(F, ZP, Z, model.H, t) # F_t = Z_t * P_t * Z_t + H
            StateSpaceModels.update_P_Ztransp_Finv!(P_Ztransp_invF, ZP, F, t) # P_Ztransp_invF   = ZP' * invF(F, t)
            StateSpaceModels.update_K!(K, P_Ztransp_invF, model.T, t) # K_t = T * P_t * Z_t * F^-1_t
            StateSpaceModels.update_att!(att, a, P_Ztransp_invF, v, t) # attt = a_t + P_t * Z_t * F^-1_t * v_t
            StateSpaceModels.update_Ptt!(Ptt, P, P_Ztransp_invF, ZP, t) # Pttt = P_t - P_t * Z_t' * F^-1_t * Z_t * P_t
            StateSpaceModels.update_a!(a, att, model.T, c, t) # a_t+1 = T * attt + c_t
            StateSpaceModels.update_P!(P, model.T, Ptt, RQR, t) # P_t+1 = T * Pttt * T' + RQR'
            if time_invariant && StateSpaceModels.check_steady_state(P, t, tol)
                steadystate = true
                tsteady     = t
            end
        end
    end
    # Return the auxiliary filter structre
    return KalmanFilter(a[model.dim.n+1:end,:], att[model.dim.n+1:end,:], v, P, Ptt, F, steadystate, tsteady, K)
end

lampspuc / statespacemodels.jl Goto Github PK