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 \(\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.12.7
Commit 6d172b025e4 (2026-08-15 08:05 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
  LLVM: libLLVM-18.1.7 (ORCJIT, skylake)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 12 virtual cores)

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.41.6"))
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

     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`

Pkg.status()

Status `~/.julia/environments/v1.12/Project.toml`
⌃ [6e4b80f9] BenchmarkTools v1.5.0
⌃ [31c24e10] Distributions v0.25.108
⌅ [682c06a0] JSON v0.21.4
⌃ [b964fa9f] LaTeXStrings v1.3.1
⌃ [d96e819e] Parameters v0.12.3
⌃ [91a5bcdd] Plots v1.41.6
⌃ [86711068] RxInfer v3.6.0
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

## 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(\lambda_x) &= \mathcal{Gam}(\lambda_x;\ a, b)\\ p(\boldsymbol{\theta}) &= \mathcal{N}(\boldsymbol{\theta};\ \mathbf{m}_{\theta}, \mathbf{V}_{\theta})\\ p(\mathbf{x}_t \mid \mathbf{x}_{t-1}, \boldsymbol{\theta}) &= \mathcal{N}\left(\mathbf{x}_t;\ \mathbf{A}\left(\boldsymbol{\theta}\right) \mathbf{x}_{t-1},\, \mathbf{V}(\lambda_x)\right)\\ p(y_t \mid \mathbf{x}_t) &= \mathcal{N}(y_t;\ \mathbf{x}_t, \lambda_y^{-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{x}_{t-1} = (x_{t-1}, x_{t-2}, ..., x_{t-M^{AR}})^\top\) is the previous state
\(\mathbf{x}_t = (x_{t}, x_{t-1}, ..., x_{t-M^{AR}+1})^\top\) is the current state
\(y_t\) is the current observation
\(\lambda_x\) is the transition precision
\(\lambda_y\) 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}(\lambda_x) = (1/ \lambda_x)\mathbf{v}_u\mathbf{v}_u^\top\) is the transition covariance
\(\mathbf{A}(\boldsymbol{\theta}) = \begin{bmatrix} \boldsymbol{\qquad\theta}^\top \\ \mathbf{I}_{M^{AR}-1} & \mathbf{0} \end{bmatrix}\) is the transition matrix

## Unicode character for ꞉ to indicate time sequences
## Copy and paste this printed form to the rest of the code identifier
## The previously used \lmrk (triangular) colon silently failed 
##  in GraphPPL.
print(Char(0xA789)) ## prints ꞉
repr(UInt32(last("_xˣ꞉"))) ## verify "0x0000a789"

꞉

"0x0000a789"

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
    # _xˣː = Vector{Vector{Float64}}(undef, N + 3Mᴬᴿ)
    _xˣ꞉ = Vector{Vector{Float64}}(undef, N + 3Mᴬᴿ)
    y꞉ = Vector{Float64}(undef, N + 3Mᴬᴿ)
    λˣₓ_std = sqrt(inv(λˣₓ)) ## transition standard dev
    λˣᵧ_std = sqrt(inv(λˣᵧ)) ## observation standard dev
    _xˣ꞉[1] = randn(rng, Mᴬᴿ) ## Mᴬᴿ random numbers

    ## +++ Matrix implementation
    ## --- 
    for t in 2:(N + 3Mᴬᴿ)
        ## TRANSITION
        ## Implementation 1: ϵ is implied in mean of Normal()
        ## _xˣ[t] = vcat(
        ##     rand( ##. noised current value is pushed in at the top of state vector
        ##         rng, 
        ##         Normal(dot(_θˣ, _xˣː[t-1]), λˣₓ_std)
        ##     ), 
        ##     _xˣː[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))
        _xˣ꞉[t] = vcat(
            dot(_θˣ, _xˣ꞉[t-1]) + ϵˣ, 
            _xˣ꞉[t-1][1:end-1] ##. rest shifted down, last element pushed out & lost
        )
        ## _xˣː[t] = rand(rng, MvNormal(BIθₜI*_xˣː[t-1], VIλₓI)) ##. Eq (2b)

        ## OBSERVATION
        y꞉[t] = rand(rng, Normal(_xˣ꞉[t][1], λˣᵧ_std))
        ## y꞉[t] = rand(rng, Normal(transpose(vᵤ)*_xˣː[t], λˣᵧ_std)) ##. Eq (2c)
    end
    
    return _xˣ꞉[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
    # _xˣː = Vector{Vector{Float64}}(undef, N + 3Mᴬᴿ)
    _xˣ꞉ = Vector{Vector{Float64}}(undef, N + 3Mᴬᴿ)
    y꞉ = Vector{Float64}(undef, N + 3Mᴬᴿ)
    λˣₓ_std = sqrt(inv(λˣₓ)) ## transition standard dev
    λˣᵧ_std = sqrt(inv(λˣᵧ)) ## observation standard dev
    _xˣ꞉[1] = randn(rng, Mᴬᴿ) ## Mᴬᴿ random numbers

    ## +++ Matrix implementation
    ## Build _AIθₜI, the transition matrix
    _Iₘ₋₁ = Diagonal(ones(Mᴬᴿ - 1))
    _V0 = Vector(zeros(Mᴬᴿ - 1)) ##vector of zeros
    _Aᵇᵒᵗᵗᵒᵐ = hcat(_Iₘ₋₁, _V0)
    _θₜ = _θˣ
    _AIθₜI = vcat(transpose(_θₜ), _Aᵇᵒᵗᵗᵒᵐ)
    ## 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()
        # _xˣː[t] = rand(rng, MvNormal(_AIθₜI*_xˣː[t-1], VIλₓI)) ##. Eq (2b)
        ## Implementation 2: ϵ is calced explicitly (aligns with definition)        
        ϵˣ = rand(rng, MvNormal(Vector(zeros(Mᴬᴿ)), _VIλₓI))
        _xˣ꞉[t] = _AIθₜI*_xˣ꞉[t-1] + ϵˣ

        ## OBSERVATION
        y꞉[t] = rand(rng, Normal(transpose(_vᵤ)*_xˣ꞉[t], λˣᵧ_std)) ##. Eq (2c)
    end
    
    return _xˣ꞉[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)

## _xˣ꞉, y꞉ = generate_sales_data_nonmatrix(rng, N, _θˣ, λˣₓ, λˣᵧ);
_xˣ꞉, y꞉ = generate_sales_data_matrix(rng, N, _θˣ, λˣₓ, λˣᵧ);

_xˣ꞉[1:7]

7-element Vector{Vector{Float64}}:
 [-53.62504471531288, -35.328336263380315]
 [-74.23708315021511, -53.62504471531507]
 [-93.27315652828227, -74.23708315021506]
 [-113.86673111848424, -93.27315652828311]
 [-135.82053836654248, -113.86673111848494]
 [-154.61663765515829, -135.82053836654092]
 [-169.84157498646235, -154.61663765515746]

y꞉[1:7]

7-element Vector{Float64}:
 -116.09952583065713
   -9.453371850211553
  -93.94444456544343
  -85.58683869381659
 -123.46833802328356
 -170.32619934639808
 -135.95594370336158

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

offset = 500.0
_xˣ꞉ = _xˣ꞉ + [[offset, offset] for i=1:N]
y꞉ = y꞉ + [offset for i=1:N];

Let’s plot our synthetic dataset:

plot(
    1:N,
    first.(_xˣ꞉) .+ offset, 
    title="Simulated sales data",
    xlabel=L"$\mathrm{time}, t$",
    color="red",
    label=L"$\mathrm{Hidden\ state\ component}, \mathbf{x}^*_t$",
) #.
        scatter!(1:N, y꞉ .+ offset, 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(
y꞉, ##. 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
        _x₀ ~ MvNormal(μ = zeros(Mᴬᴿ), Λ = diageye(Mᴬᴿ)) ##.kw μ,Λ only work inside macro
    else ## Univariate
        _θ  ~ Normal(μ = 0.0, γ = 1.0) ##. γ is RxInfer's precision keyword, not the r.v. λₓ
        _x₀ ~ Normal(μ = 0.0, γ = 1.0)
    end
    _xₜ₋₁ = _x₀
    for t in eachindex(y꞉)
        _x꞉[t] ~ AR(x = _xₜ₋₁, θ = _θ, γ = λₓ) #.Eq (2b); x not _xː (GraphPPL)
        if 𝚃ᴬᴿ === Multivariate
            y꞉[t] ~ Normal(μ = dot(_vᵤ, _x꞉[t]), γ = λᵧ) #.Eq (2c)
        else
            y꞉[t] ~ Normal(μ = _vᵤ*_x꞉[t], γ = λᵧ) #.Eq (2c)
        end
        _xₜ₋₁ = _x꞉[t]
    end
end

@constraints function ar_constraints() 
    q(_x₀, _x꞉, _θ, λₓ) = q(_x₀, _x꞉)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 = (y꞉=y꞉,), #.
            constraints   = ar_constraints(),
            meta          = ar_meta(𝚃ᴬᴿ, m, ARsafe()),
            options       = (limit_stack_depth = 100,),
            initialization = ar_init(𝚃ᴬᴿ, m),
            returnvars    = (_x꞉=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 (λₓ, _x꞉, _θ)
  Free Energy:     | Real[44997.2, 44981.3, 44931.5, 44593.6, 29002.6, 4928.99, 2209.64, 1732.64, 1667.44, 1653.82  …  1649.62, 1649.62, 1649.62, 1649.62, 1649.62, 1649.62, 1649.62, 1649.62, 1649.62, 1649.62]

 Inference results:
  Posteriors       | available for (λₓ, _x꞉, _θ)
  Free Energy:     | Real[45010.2, 44996.0, 44950.4, 44531.8, 17244.3, 2722.87, 1978.89, 1867.18, 1803.75, 1764.66  …  1652.54, 1652.54, 1652.54, 1652.54, 1652.54, 1652.54, 1652.54, 1652.54, 1652.54, 1652.54]

 Inference results:
  Posteriors       | available for (λₓ, _x꞉, _θ)
  Free Energy:     | Real[45020.1, 45005.9, 44966.3, 44553.6, 17672.5, 3157.61, 1947.23, 1827.08, 1811.25, 1802.91  …  1646.21, 1646.2, 1646.19, 1646.18, 1646.17, 1645.82, 1645.5, 1645.47, 1645.9, 1646.32]

 Inference results:
  Posteriors       | available for (λₓ, _x꞉, _θ)
  Free Energy:     | Real[45029.2, 45013.5, 44979.2, 44597.6, 20524.2, 4140.52, 1981.19, 1887.73, 1870.58, 1859.91  …  1647.97, 1650.4, 1651.93, 1643.26, 1636.49, 1638.42, 1629.09, 1638.96, 1642.45, 1644.45]

 Inference results:
  Posteriors       | available for (λₓ, _x꞉, _θ)
  Free Energy:     | Real[45038.2, 45019.9, 44989.8, 44645.2, 23977.0, 5100.34, 2104.49, 1958.58, 1930.85, 1901.38  …  1643.69, 1646.7, 1639.02, 1633.06, 1650.07, 1641.31, 1653.34, 1641.72, 1630.95, 1648.04]

 Inference results:
  Posteriors       | available for (λₓ, _x꞉, _θ)
  Free Energy:     | Real[45047.4, 45025.6, 44998.8, 44690.5, 27427.9, 5923.78, 2351.46, 1991.97, 1969.83, 1925.64  …  1637.57, 1642.86, 1627.7, 1661.37, 1669.77, 1619.89, 1636.73, 1618.56, 1648.52, 1652.11]

 Inference results:
  Posteriors       | available for (λₓ, _x꞉, _θ)
  Free Energy:     | Real[45056.7, 45031.0, 45006.8, 44731.9, 30507.0, 6641.37, 2822.4, 2042.54, 1984.49, 1951.36  …  1691.12, 1670.09, 1663.56, 1662.9, 1672.69, 1650.56, 1674.16, 1658.73, 1650.11, 1662.46]

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}:
 1649.6247542602655
 1652.5413450037695
 1646.3178693260152
 1644.4520667386723
 1648.0433087113252
 1652.1076666725967
 1662.46321011273

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",
        ## ylim=(0, 7)
    )
    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 $\lambda_x^*$"
    )
    return p
end

plot_transition_precision (generic function with 1 method)

plot_transition_precision(results)

# Mᵇᵉˢᵗ = 5 ## best order
Mᵇᵉˢᵗ = 4 ## best order
@unpack _x꞉, _θ, λₓ = results[Mᵇᵉˢᵗ].posteriors

Dict{Symbol, Vector} with 3 entries:
  :λₓ  => GammaShapeRate{Float64}[GammaShapeRate{Float64}(a=151.0, b=35.3412), …
  :_x꞉ => MvNormalWeightedMeanPrecision{Float64, Vector{Float64}, Matrix{Float6…
  :_θ  => MvNormalWeightedMeanPrecision{Float64, Vector{Float64}, Matrix{Float6…

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

p2 = plot(
    mean.(λₓ), ribbon=std.(λₓ), 
    fillalpha=0.2,
    label=L"Inferred transition precision $\lambda_x$", 
    color="green",
    legend=:topright,
    xlabel="iterations",
)
p2 = plot!(
    [λˣₓ], 
    seriestype=:hline, 
    color="red", #.
    label=L"Real transition precision $\lambda_x^*$"
)

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.(_xˣ꞉)[subrange], 
    title=L"Subrange of: AR, $M^{AR}=$"*"$Mᵇᵉˢᵗ",
    xlabel=L"time $t$",
    color="red",
    label=L"Hidden state component $\mathbf{x}^*_t$")
scatter!(subrange, y꞉[subrange], label=L"Observation $x_t=y_t$", color="orange")
plot!(
    subrange, 
    first.(mean.(_x꞉))[subrange],
    ribbon=sqrt.(first.(var.(_x꞉)))[subrange],
    fillalpha=0.2,
    label=L"Inferred state component $\mathbf{x}_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])

4-element Vector{Float64}:
  1.045257779867135
  0.08830020968205712
  0.09032250642577659
 -0.2457606622600396

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

(1.045257779867135, 0.08830020968205712, 0.09032250642577659, -0.2457606622600396)

We decide to forecast 10 values into the future by having \(n^{\mathrm{fcast}}=10\).

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

x꞉ᶠᶜᵃˢᵗ = Vector{Float64}(undef, nᶠᶜᵃˢᵗ);

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

310-element Vector{Float64}:
 272.94080430289796
 389.0984769577544
 390.49634157427784
 415.42391521448025
 389.2852440214288
 356.40229429790577
 357.6686366712435
 334.8959921546506
 324.43930927760925
 304.58982534026444
   ⋮
   2.060710655e-315
   9.34e-322
   5.0e-324
   6.94607119497607e-310
   9.4e-322
   6.94607119497765e-310
   6.9461272312858e-310
   9.44e-322
   9.999999999999965e-13

Now we calculate the forecast values:

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

310-element Vector{Float64}:
 272.94080430289796
 389.0984769577544
 390.49634157427784
 415.42391521448025
 389.2852440214288
 356.40229429790577
 357.6686366712435
 334.8959921546506
 324.43930927760925
 304.58982534026444
   ⋮
 722.2133846233114
 704.9939444277667
 685.8196845599591
 664.5875387649759
 641.4084976251897
 617.8056154343553
 593.8823431826843
 569.9166544146792
 546.318526846662

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