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

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 \(\boldsymbol\theta\) of the AR process. The inference process also learned the value of the transition precision parameter \(\gamma\) 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(\gamma) &= \mathcal{Gam}(\gamma|a, b)\\ p(\boldsymbol{\theta}) &= \mathcal{N}(\boldsymbol{\theta}|\mathbf{m}_{\theta}, \mathbf{V}_{\theta})\\ p(\mathbf{s}_t \mid \mathbf{s}_{t-1}, \boldsymbol{\theta}) &= \mathcal{N}\left(\mathbf{s}_t \mid \mathbf{B}\left(\boldsymbol{\theta}\right) \mathbf{s}_{t-1},\, \mathbf{V}(\gamma)\right)\\ p(x_t \mid \mathbf{s}_t) &= \mathcal{N}(x_t \mid \mathbf{s}_t, \tau^{-1}) \end{aligned}\]

where

\(M^{AR}\) is the order of the autoregression process
\(\boldsymbol{\theta}=(\theta_1, \theta_2, ..., \theta_{M^{AR}})^{\top}\) is a vector of transition coefficients
\(\mathbf{s}_{t-1} = (s_{t-1}, s_{t-2}, ..., s_{t-M^{AR}})^\top\) is the previous state
\(\mathbf{s}_t = (s_{t}, s_{t-1}, ..., s_{t-M^{AR}+1})^\top\) is the current state
\(\gamma\) is the transition precision
\(\tau\) is the observation precision
\(\mathbf{v}_u = (1,0,...,0)^\top\) is an \(M^{AR}\)-dimensional unit vector selecting the first component of the vector that follows
\(\mathbf{V}(\gamma) = (1/ \gamma)\mathbf{v}_u\mathbf{v}_u^\top\) is the transition covariance
\(\mathbf{B}(\boldsymbol{\theta}) = \begin{bmatrix} \boldsymbol{\qquad\theta}^\top \\ \mathbf{I}_{M^{AR}-1} & \mathbf{0} \end{bmatrix}\) is the transition matrix

To generate some synthetic sales data, we use a second-order autoregressive (AR) process, i.e. \(M^{AR}=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^{AR} > 1\)
Univariate AR (reduces to a simple State-Space-Model), with order \(M^{AR}\) = 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^{AR}=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 \(\theta\) and call it \(\theta^{\mathrm{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^{\mathrm{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$")