Sales Forecasting with Time-Varying Autoregressive Models (Part 3)

Back to Portfolio of Projects | LearnableLoopAI.com | Blog |

In Part 2 of the Sales Forecasting project we setup the foundations of the design and provided a simple implementation. In Part 3 we do the first forecast beyond the available synthetic data. As before, we used RxInfer to learn the best order (5) as well as the parameters $θ$ of the AR process. The inference process also learned the value of the transition precision parameter $γ$ as well as the $s_{t}$ state component.

versioninfo() ## Julia version

Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × Intel(R) Core(TM) i7-8700B CPU @ 3.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 12 virtual cores)
Environment:
  JULIA_NUM_THREADS =

import Pkg
Pkg.add(Pkg.PackageSpec(;name="RxInfer", version="3.6.0"))
Pkg.add(Pkg.PackageSpec(;name="Distributions", version="0.25.108"))
Pkg.add(Pkg.PackageSpec(;name="Random"))
Pkg.add(Pkg.PackageSpec(;name="Plots", version="1.40.1"))
Pkg.add(Pkg.PackageSpec(;name="LaTeXStrings", version="1.3.1"))
Pkg.add(Pkg.PackageSpec(;name="BenchmarkTools", version="1.5.0"))
Pkg.add(Pkg.PackageSpec(;name="Parameters", version="0.12.3"))

using RxInfer, Distributions, LinearAlgebra, Random, Plots, LaTeXStrings, BenchmarkTools, Parameters

   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`
   Resolving package versions...
  No Changes to `~/.julia/environments/v1.10/Project.toml`
  No Changes to `~/.julia/environments/v1.10/Manifest.toml`

Pkg.status()

Status `~/.julia/environments/v1.10/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
⌃ [31c24e10] Distributions v0.25.108
  [b964fa9f] LaTeXStrings v1.3.1
  [d96e819e] Parameters v0.12.3
⌃ [91a5bcdd] Plots v1.40.1
  [86711068] RxInfer v3.6.0
  [9a3f8284] Random
Info Packages marked with ⌃ have new versions available and may be upgradable.

## macro h(x) ## for help
#     quote
#         display("text/markdown", @doc $x)
#     end    
# end
# @h Normal() ## help on Normal()

Our model for the Variational Bayesian Inference of a latent autoregressive model can be represented mathematically as follows:

$\begin{aligned} p (γ) & = Gam (γ | a, b) \\ p (θ) & = N (θ | m_{θ}, V_{θ}) \\ p (s_{t} ∣ s_{t - 1}, θ) & = N (s_{t} ∣ B (θ) s_{t - 1}, V (γ)) \\ p (x_{t} ∣ s_{t}) & = N (x_{t} ∣ s_{t}, τ^{- 1}) \end{aligned}$

where

$M^{A R}$ is the order of the autoregression process
$θ = (θ_{1}, θ_{2}, . . ., θ_{M^{A R}})^{⊤}$ is a vector of transition coefficients
$s_{t - 1} = (s_{t - 1}, s_{t - 2}, . . ., s_{t - M^{A R}})^{⊤}$ is the previous state
$s_{t} = (s_{t}, s_{t - 1}, . . ., s_{t - M^{A R} + 1})^{⊤}$ is the current state
$γ$ is the transition precision
$τ$ is the observation precision
$v_{u} = (1, 0, . . ., 0)^{⊤}$ is an $M^{A R}$ -dimensional unit vector selecting the first component of the vector that follows
$V (γ) = (1 / γ) v_{u} v_{u}^{⊤}$ is the transition covariance
$B (θ) = [\begin{matrix} θ^{⊤} \\ I_{M^{A R} - 1} & 0 \end{matrix}]$ is the transition matrix

To generate some synthetic sales data, we use a second-order autoregressive (AR) process, i.e. $M^{A R} = 2$ with the following parameters:

_f = .02 ## Frequency of a sine wave (Hz)
_ω = 2π*_f
_θˢⁱⁿᵉ = [2cos(_ω), -1]

2-element Vector{Float64}:
  1.9842294026289558
 -1.0

function generate_sales_data_nonmatrix(
rng, ## random number generator
N, ## number of data points 
θ̃, ## AR parameters
γ̃, ## transition precision
τ̃) ## observation precision
    Mᴬᴿ = length(θ̃) ## order of process
    s̃ = Vector{Vector{Float64}}(undef, N + 3Mᴬᴿ)
    y = Vector{Float64}(undef, N + 3Mᴬᴿ)
    γ̃_std = sqrt(inv(γ̃)) ## transition standard dev
    τ̃_std = sqrt(inv(τ̃)) ## observation standard dev
    s̃[1] = randn(rng, Mᴬᴿ) ## Mᴬᴿ random numbers

    ## +++ Matrix implementation
    ## --- 
    for t in 2:(N + 3Mᴬᴿ)
        ## TRANSITION
        ## Implementation 1: ϵ is implied in mean of Normal()
        ## s̃[t] = vcat(
        ##     rand( ##. noised current value is pushed in at the top of state vector
        ##         rng, 
        ##         Normal(dot(θ̃, s̃[t-1]), γ̃_std)
        ##     ), 
        ##     s̃[t-1][1:end-1] ##. rest shifted down, last element pushed out & lost
        ## )
        ## Implementation 2: ϵ is calced explicitly (aligns with definition)
        ϵ̃ = rand(rng, Normal(0, γ̃_std))
        s̃[t] = vcat(
            dot(θ̃, s̃[t-1]) + ϵ̃, 
            s̃[t-1][1:end-1] ##. rest shifted down, last element pushed out & lost
        )
        ## s̃[t] = rand(rng, MvNormal(BIθₜI*s̃[t-1], VIγI)) ##. Eq (2b)

        ## OBSERVATION
        y[t] = rand(rng, Normal(s̃[t][1], τ̃_std))
        ## y[t] = rand(rng, Normal(transpose(vᵤ)*s̃[t], τ̃_std)) ##. Eq (2c)
    end
    
    return s̃[1+3Mᴬᴿ:end], y[1+3Mᴬᴿ:end] #.
end

generate_sales_data_nonmatrix (generic function with 1 method)

function generate_sales_data_matrix(
rng, ## random number generator
N, ## number of data points 
θ̃, ## AR parameters
γ̃, ## transition precision
τ̃) ## observation precision
    Mᴬᴿ = length(θ̃) ## order of process
    s̃ = Vector{Vector{Float64}}(undef, N + 3Mᴬᴿ)
    y = Vector{Float64}(undef, N + 3Mᴬᴿ)
    γ̃_std = sqrt(inv(γ̃)) ## transition standard dev
    τ̃_std = sqrt(inv(τ̃)) ## observation standard dev
    s̃[1] = randn(rng, Mᴬᴿ) ## Mᴬᴿ random numbers

    ## +++ Matrix implementation
    ## Build BIθₜI, the transition matrix
    Iₘ₋₁ = Diagonal(ones(Mᴬᴿ - 1))
    V0 = Vector(zeros(Mᴬᴿ - 1)) ##vector of zeros
    Bᵇᵒᵗᵗᵒᵐ = hcat(Iₘ₋₁, V0)
    θₜ = θ̃
    BIθₜI = vcat(transpose(θₜ), Bᵇᵒᵗᵗᵒᵐ)
    ## Build VIγI, the transition covariance matrix
    vᵤ = Vector(vcat(1.0, zeros(Mᴬᴿ - 1)))
    VIγI = (1/γ̃)*vᵤ*transpose(vᵤ)
    for i in 2:Mᴬᴿ
        VIγI[i, i] = 1e-12*tiny ##tiny values on remaining off-diagonals to avoid PosDefException
    end
    ## --- 
    for t in 2:(N + 3Mᴬᴿ)
        ## TRANSITION
        ## Implementation 1: ϵ is implied in mean of Normal()
        # s̃[t] = rand(rng, MvNormal(BIθₜI*s̃[t-1], VIγI)) ##. Eq (2b)
        ## Implementation 2: ϵ is calced explicitly (aligns with definition)        
        ϵ̃ = rand(rng, MvNormal(Vector(zeros(Mᴬᴿ)), VIγI))
        s̃[t] = BIθₜI*s̃[t-1] + ϵ̃

        ## OBSERVATION
        y[t] = rand(rng, Normal(transpose(vᵤ)*s̃[t], τ̃_std)) ##. Eq (2c)
    end
    
    return s̃[1+3Mᴬᴿ:end], y[1+3Mᴬᴿ:end] #.
end

generate_sales_data_matrix (generic function with 1 method)

## Seed for reproducibility
_seed = 123
_rng  = MersenneTwister(_seed)

## Number of observations in synthetic dataset
_N = 300

## AR process parameters
_γ̃ = 0.1 ## real γ (transition precision)

_τ̃ = 0.001 ## real τ (observation precision)
_θ̃ = _θˢⁱⁿᵉ ## real θ (AR parameters)

## _s̃s, _ys = generate_sales_data_nonmatrix(_rng, _N, _θ̃, _γ̃, _τ̃);
_s̃s, _ys = generate_sales_data_matrix(_rng, _N, _θ̃, _γ̃, _τ̃);

First, we offset the generated values to make the sales data positive:

offset = 500.0
_s̃s = _s̃s + [[offset, offset] for i=1:_N]
_ys = _ys + [offset for i=1:_N];

Let’s plot our synthetic dataset:

plot(
    1:_N,
    first.(_s̃s), 
    title="Simulated sales data",
    xlabel=L"$\mathrm{time}, t$",
    color="red",
    label=L"$\mathrm{Hidden\ state\ component}, \tilde{s}_t$",
) #.
scatter!(1:_N, _ys, label=L"$\mathrm{Observation}, y_t$", color="orange")

Now we specify a probabilistic model, the inference constraints and then run an inference procedure with RxInfer. We need two different models for:

Multivariate AR, with order $M^{A R} > 1$
Univariate AR (reduces to a simple State-Space-Model), with order $M^{A R}$ = 1

@model function lar_model(
x, ##. data/observations 
𝚃ᴬᴿ, ##. Uni/Multi variate 
Mᴬᴿ, ##. AR order
vᵤ, ##. unit vector 
τ) ##. observation precision     
    ## Priors
    γ  ~ Gamma(α = 1.0, β = 1.0) ##. for transition precision    
    if 𝚃ᴬᴿ === Multivariate
        θ  ~ MvNormal(μ = zeros(Mᴬᴿ), Λ = diageye(Mᴬᴿ)) ##.kw μ,Λ only work inside macro
        s₀ ~ MvNormal(μ = zeros(Mᴬᴿ), Λ = diageye(Mᴬᴿ)) ##.kw μ,Λ only work inside macro
    else ## Univariate
        θ  ~ Normal(μ = 0.0, γ = 1.0)
        s₀ ~ Normal(μ = 0.0, γ = 1.0)
    end
    sₜ₋₁ = s₀
    for t in eachindex(x)
        s[t] ~ AR(sₜ₋₁, θ, γ) #.Eq (2b)
        if 𝚃ᴬᴿ === Multivariate
            x[t] ~ Normal(μ = dot(vᵤ, s[t]), γ = τ) #.Eq (2c)
        else
            x[t] ~ Normal(μ = vᵤ*s[t], γ = τ) #.Eq (2c)
        end
        sₜ₋₁ = s[t]
    end
end

@constraints function ar_constraints() 
    q(s₀, s, θ, γ) = q(s₀, s)q(θ)q(γ)
end

ar_constraints (generic function with 1 method)

@meta function ar_meta(𝚃ᴬᴿ, Mᴬᴿ, stype)
    AR() -> ARMeta(𝚃ᴬᴿ, Mᴬᴿ, stype)
end

ar_meta (generic function with 1 method)

@initialization function ar_init(𝚃ᴬᴿ, Mᴬᴿ)
    q(γ) = GammaShapeRate(1.0, 1.0)
    if 𝚃ᴬᴿ === Multivariate    
        q(θ) = MvNormalMeanPrecision(zeros(Mᴬᴿ), diageye(Mᴬᴿ))
    else ## Univariate
        q(θ) = NormalMeanPrecision(0.0, 1.0)
    end
end

ar_init (generic function with 1 method)

function run_inference(Nmodels)
    results = Vector{Any}(undef, Nmodels)
    for m = 1:Nmodels
        println("m = $m")
        if m > 1 𝚃ᴬᴿ = Multivariate elseif m === 1 𝚃ᴬᴿ = Univariate
        else error("m=$m which is not an allowed value.") end    
        results[m] = infer(
            model = lar_model(
                𝚃ᴬᴿ=𝚃ᴬᴿ, ##. Uni/Multi variate
                Mᴬᴿ=m, ##. AR order
                vᵤ=ReactiveMP.ar_unit(𝚃ᴬᴿ, m), ##. unit vector
                τ=_τ̃), ##. observation precision
            data = (x=_ys,), #.
            constraints   = ar_constraints(),
            meta          = ar_meta(𝚃ᴬᴿ, m, ARsafe()),
            options       = (limit_stack_depth = 100,),
            initialization = ar_init(𝚃ᴬᴿ, m),
            returnvars    = (s=KeepLast(), γ=KeepEach(), θ=KeepEach()),
            free_energy   = true,
            iterations    = 100,
            showprogress  = false
        )
    end
    return(results)
end

run_inference (generic function with 1 method)

Now we run the inference for $M^{A R} = 1, 2, . . ., 7$ :

## _results = run_inference(length(_θ̃))
_results = run_inference(length(_θ̃) + 5)

m = 1
m = 2
m = 3
m = 4
m = 5
m = 6
m = 7

7-element Vector{Any}:
 Inference results:
  Posteriors       | available for (γ, s, θ)
  Free Energy:     | Real[42701.1, 42687.1, 42645.9, 42398.3, 33951.0, 5623.04, 2468.23, 1780.97, 1686.3, 1666.36  …  1660.06, 1660.06, 1660.06, 1660.06, 1660.06, 1660.06, 1660.06, 1660.06, 1660.06, 1660.06]

 Inference results:
  Posteriors       | available for (γ, s, θ)
  Free Energy:     | Real[42713.4, 42701.1, 42664.0, 42371.7, 24685.8, 3932.85, 2154.2, 1980.86, 1894.82, 1835.54  …  1666.49, 1666.49, 1666.49, 1666.49, 1666.49, 1666.49, 1666.49, 1666.49, 1666.49, 1666.49]

 Inference results:
  Posteriors       | available for (γ, s, θ)
  Free Energy:     | Real[42722.9, 42710.6, 42679.0, 42403.7, 24323.7, 3926.1, 2166.01, 1921.63, 1890.66, 1885.2  …  1661.8, 1661.8, 1661.8, 1661.8, 1661.8, 1661.8, 1661.8, 1661.8, 1661.8, 1661.8]

 Inference results:
  Posteriors       | available for (γ, s, θ)
  Free Energy:     | Real[42731.7, 42718.0, 42691.1, 42448.6, 26840.2, 4854.59, 2293.54, 2010.31, 1984.34, 1960.02  …  1657.13, 1657.15, 1657.15, 1657.16, 1657.17, 1657.17, 1657.17, 1657.18, 1657.18, 1657.18]

 Inference results:
  Posteriors       | available for (γ, s, θ)
  Free Energy:     | Real[42740.5, 42724.3, 42701.1, 42492.4, 29895.0, 5858.92, 2753.88, 2099.35, 2055.19, 2004.08  …  1669.54, 1666.56, 1640.63, 1643.82, 1664.04, 1666.34, 1670.33, 1662.44, 1666.39, 1639.2]

 Inference results:
  Posteriors       | available for (γ, s, θ)
  Free Energy:     | Real[42749.4, 42729.9, 42709.7, 42531.5, 32672.1, 6762.14, 3615.16, 2182.45, 2116.36, 2058.74  …  1673.94, 1699.09, 1731.72, 1708.24, 1688.49, 1712.02, 1701.23, 1687.45, 1667.2, 1688.61]

 Inference results:
  Posteriors       | available for (γ, s, θ)
  Free Energy:     | Real[42758.5, 42735.2, 42717.3, 42565.2, 34951.9, 7668.57, 4652.86, 2307.96, 2166.63, 2096.53  …  1711.05, 1716.59, 1696.4, 1698.84, 1726.57, 1662.98, 1658.88, 1671.77, 1683.77, 1670.59]

function plot_free_energy(results)
    p = plot(title="Bethe Free Energy", legend=:topright, xlabel="iterations")
    for m = 1:length(results)
        p = plot!(p, results[m].free_energy, label="AR$m")
    end
    return p
end

plot_free_energy (generic function with 1 method)

plot_free_energy(_results)

values = [_results[m].free_energy[end] for m in 1:length(_results)]

7-element Vector{Float64}:
 1660.062561707014
 1666.4864438734899
 1661.7964138139591
 1657.1795897457341
 1639.1966689677083
 1688.6112655905781
 1670.590204262484

function plot_converged_free_energy(results)
    values = [results[m].free_energy[end] for m in 1:length(results)]
    bar(
        1:length(results), 
        values, 
        legend=:none, 
        ylim=(1600, 1700), 
        title="Converged Free Energy vs AR Order"
    )
    xlabel!("AR order")
    ylabel!("Minimum Free Energy")
end

plot_converged_free_energy (generic function with 1 method)

plot_converged_free_energy(_results)

function plot_transition_precision(results)
    p = plot(title="Inferred transition precision", legend=:topright, xlabel="iterations")
    for m = 1:length(results)
        p = plot!(p, mean.(results[m].posteriors[:γ]), label="AR$m")
    end
    plot!(
        p, [_γ̃], seriestype=:hline, color="red", linestyle=:dash, 
        label=L"Real transition precision $\tilde{γ}$"
    )
    return p
end

plot_transition_precision (generic function with 1 method)

plot_transition_precision(_results)

Mᵇᵉˢᵗ = 5 ## best order
@unpack s, θ, γ = _results[Mᵇᵉˢᵗ].posteriors

Dict{Symbol, Vector} with 3 entries:
  :γ => GammaShapeRate{Float64}[GammaShapeRate{Float64}(a=151.0, b=30.1702), Ga…
  :s => MvNormalWeightedMeanPrecision{Float64, Vector{Float64}, Matrix{Float64}…
  :θ => MvNormalWeightedMeanPrecision{Float64, Vector{Float64}, Matrix{Float64}…

_p1 = plot(
    first.(_s̃s), ## first.() gives s̃ₜ, not the delayed s̃ₜ₋₁, s̃ₜ₋₂, ..., s̃ₜ₋ₘ₊₁
    title="AR, " * L"$M^{AR}=$" * "$Mᵇᵉˢᵗ",
    xlabel=L"time $t$",
    color="red",
    label=L"Hidden state component $\tilde{s}_t$")
_p1 = scatter!(_p1, _ys, label=L"Observation $x_t=y_t$", color="orange")
_p1 = plot!(
    _p1, first.(mean.(s)), 
    ribbon=first.(std.(s)), 
    fillalpha=0.2,
    label=L"Inferred state component $s_t$",
    color="green",
    legend=:bottomleft)

_p2 = plot(
    mean.(γ), ribbon=std.(γ), 
    fillalpha=0.2,
    label=L"Inferred transition precision $\gamma$", 
    color="green",
    legend=:topright,
    xlabel="iterations",
)
_p2 = plot!(
    [_γ̃], 
    seriestype=:hline, 
    color="red", #.
    label=L"Real transition precision $\tilde{γ}$"
)

_p3 = plot(_results[Mᵇᵉˢᵗ].free_energy, label="Bethe Free Energy", xlabel="iterations")
plot(_p1, _p2, _p3, layout=@layout([ a; b c ]))

Let’s look at a smaller range of data for the sake of interest:

_subrange = 200 : _N
plot(
    _subrange, 
    first.(_s̃s)[_subrange], 
    title=L"Subrange of: AR, $M^{AR}=$"*"$Mᵇᵉˢᵗ",
    xlabel=L"time $t$",
    color="red",
    label=L"Hidden state component $\tilde{s}_t$")
scatter!(_subrange, _ys[_subrange], label=L"Observation $x_t=y_t$", color="orange")
plot!(
    _subrange, 
    first.(mean.(s))[_subrange],
    ribbon=sqrt.(first.(var.(s)))[_subrange],
    fillalpha=0.2,
    label=L"Inferred state component $s_t$",
    color="green",
    legend=:bottomleft)

Finally, we look at how the learned AR coefficients converge:

let
    pθ = plot()
    θms = mean.(θ) ##. θ means
    θvs = var.(θ) ##. θ variances
    l = length(θms)
    edim(e) = (a) -> map(r -> r[e], a)
    for t in 1:length(first(θms)) #.
        tmp = latexstring("\$ \\theta_{$(t),t} \$")
        pθ = plot!(
            pθ, 
            θms |> edim(t),
            ribbon=θvs |> edim(t) .|> sqrt, 
            fillalpha=0.2,
            label="Inferred $tmp")
    end
    for t in 1:length(_θ̃)
        tmp = latexstring("\$ \\theta_{$(t),t} \$")
        pθ = plot!(
            pθ, 
            [ _θ̃[t] ],
            seriestype=:hline, 
            label="Real $tmp")
    end
    plot(
        pθ, 
        title=L"AR coefficients, $M^{AR}=$"*"$Mᵇᵉˢᵗ",
        xlabel=L"time $t$",
        legend=:outertopright, 
        size=(800, 450)) #.
end

Finally, we forecast a number of values into the future. First, we find the best value for $θ$ and call it $θ^{best}$ .

_θᵇᵉˢᵗ = mean(θ[end])

5-element Vector{Float64}:
  0.9444189816485971
  0.20485343906764922
 -0.1213157360930712
  0.10833466838389029
 -0.18507076644769016

_θᵇᵉˢᵗ[1], _θᵇᵉˢᵗ[2], _θᵇᵉˢᵗ[3], _θᵇᵉˢᵗ[4], _θᵇᵉˢᵗ[5]

(0.9444189816485971, 0.20485343906764922, -0.1213157360930712, 0.10833466838389029, -0.18507076644769016)

We decide to forecast 10 values into the future by having $n^{fcast} = 10$ .

_nᶠᶜᵃˢᵗ = 10 ## number of forecast values

_sᶠᶜᵃˢᵗ = Vector{Float64}(undef, _nᶠᶜᵃˢᵗ);

## concat existing data with forecast to get complete signal
_sᶜᵒᵐᵖˡ = [first.(mean.(s)); _sᶠᶜᵃˢᵗ]

310-element Vector{Float64}:
 341.3856591972407
 428.7163763705686
 499.00226325690613
 483.8842281777904
 492.14833325413315
 453.1403215504383
 445.77905407498463
 438.9422640092849
 466.0210265298952
 443.83927752426086
   ⋮
   6.9027617478467e-310
   6.90276594706e-310
   6.9027650371239e-310
   6.9026213185935e-310
   6.9027617478467e-310
   6.90276594706e-310
   6.9027650371239e-310
   6.90276147999355e-310
   6.9027617478467e-310

Now we calculate the forecast values:

for t in _N+1:_N+_nᶠᶜᵃˢᵗ
    _sᶜᵒᵐᵖˡ[t] = _θᵇᵉˢᵗ[1]*_sᶜᵒᵐᵖˡ[t-1] + _θᵇᵉˢᵗ[2]*_sᶜᵒᵐᵖˡ[t-2] + _θᵇᵉˢᵗ[3]*_sᶜᵒᵐᵖˡ[t-3] + _θᵇᵉˢᵗ[4]*_sᶜᵒᵐᵖˡ[t-4]
end
_sᶜᵒᵐᵖˡ

310-element Vector{Float64}:
 341.3856591972407
 428.7163763705686
 499.00226325690613
 483.8842281777904
 492.14833325413315
 453.1403215504383
 445.77905407498463
 438.9422640092849
 466.0210265298952
 443.83927752426086
   ⋮
 362.8198419583558
 400.05132270786595
 445.00652349546334
 493.3888121836583
 547.9003228772283
 607.872993367997
 674.679865997366
 748.6875275685808
 830.8972433233798

## _subrange = 1:_N
_subrange = 200:_N
_pc1 = plot(
    _subrange,
    first.(_s̃s)[_subrange], ## first.() gives s̃ₜ, not the delayed s̃ₜ₋₁, s̃ₜ₋₂, ..., s̃ₜ₋ₘ₊₁
    title="AR, " * L"$M^{AR}=$" * "$Mᵇᵉˢᵗ",
    xlabel=L"time $t$",
    color="red",
    label=L"Hidden state component $\tilde{s}_t$",
    legend=:bottomleft)
_pc1 = scatter!(_pc1, 
    _subrange,
    _ys[_subrange], label=L"Observation $x_t=y_t$", color="orange")
_pc1 = plot!(_pc1, 
    _subrange,
    first.(mean.(s))[_subrange], 
    ribbon=first.(std.(s))[_subrange], 
    fillalpha=0.2,
    label=L"Inferred state component $s_t$",
    color="green")
_pc1 = plot!(_pc1, 
    _N+1:_N+_nᶠᶜᵃˢᵗ,
    _sᶜᵒᵐᵖˡ[_N+1:_N+_nᶠᶜᵃˢᵗ],
    title="AR, " * L"$M^{AR}=$" * "$Mᵇᵉˢᵗ",
    xlabel=L"time $t$",
    color="green",
    linestyle=:dash,
    label=L"Forecast state component $s^{\mathrm{fcast}}_t$")