How to Incorporate Explanatory Variables in State Equations Using the bnlssm Package?

Hello everyone,

I’m currently working with the bnlssm package in MATLAB to model a nonlinear state-space system. My goal is to include explanatory variables in the state equation, but I’m encountering some challenges and need your help.

Specifically, I want to modify the state equation to include external explanatory variables , such as:

And the observation equation is

In this setup: is the state variable. is an explanatory variable (could be a vector or scalar).

I have been working on implementing a time-varying A matrix to address this problem, but unfortunately, I haven’t been able to achieve success. Below, I am sharing my current code and data in the hope that you might help me identify the issue and provide a solution. Thank you so much!

data = readtable('/Users/xiaoxuan/Desktop/Matlab Code/SSE1.csv');
head(data);

SSE = data.close;
dts = data.time;
dts = datetime(dts, 'InputFormat','MM/dd/yyyy');
T = numel(SSE);
T1 = T-1;
retsp500 = price2ret(SSE);
y = retsp500 - mean(retsp500);
retdts = dts(2:end);

Z = data.high(2:end);


PriorMdl = bnlssm(@(theta)paramMap(theta, y, T1, Z), @flatLogPrior, ObservationForm="distribution", ...
    Multipoint=["A" "LogY"]);
theta0 = [0.2 0 0.5 0.7 0 1]; % Adjusted for single predictor
lower = [-1; -Inf; 0; 0; -Inf; -Inf];
upper = [1; Inf; Inf; Inf; Inf; Inf];

burnin = 1e4;
thin = 25;
rng(1)

PosteriorMdl = estimate(PriorMdl, y, theta0, Lower=lower, Upper=upper, ...
    NumParticles=500, Hessian="opg", SortParticles=false, BurnIn=burnin, Thin=thin);

function [A, B, LogY, Mean0, Cov0, StateType] = paramMap(theta,T1,Z)

    A = cell(T1, 1);  
    for t = 1:T1-1
        A{t} = @(x) theta(1) + theta(2) .* x + ...
                           theta(3) * theta(4) * exp(-0.5 .* x) .* (Z(t)-theta(5));
    end
    A{T1} = @(x) theta(1) + theta(2) .* x + ...
                   theta(3) * theta(4) * exp(-0.5 .* x) .* (Z(T1) - theta(5));

    B = theta(4) * sqrt(1 - theta(3)^2);

    LogY = @(y, x) -0.5 .* log(2*pi) - 0.5 .* x - ...
                   0.5 .* ((y - theta(5)).^2) ./ exp(x);

    Mean0 = theta(2) / (1 - theta(1));          
    Cov0 = (theta(4)^2) / (1 - theta(1)^2);

    StateType = 0;  
end
    

function logprior = flatLogPrior(theta)
    paramcon = zeros(numel(theta), 1);
    paramcon(1) = abs(theta(1)) >= 1 - eps;
    paramcon(2) = ~isfinite(theta(2)); 
    paramcon(3) = abs(theta(3)) > 1;
    paramcon(4) = theta(4) <= eps;
    paramcon(5) = ~isfinite(theta(5)); 
    
    if sum(paramcon) > 0
        logprior = -Inf;
    else
        logprior = 0;
    end
end