Does Bayesian Non-linear SSM (bnlssm) actually support time-varying coefficients?

Accodring to its documentation, 'bnlssm' object allows any of the coefficient matrices to be time-varying and, in that case, defined as cell arrays.

However, when i modify a working example given in the documentation to feature cell-arrays for either of the matrices, i get an error message that eminates from the validation check within the 'bnlssm' object (lines 700--707) failing to recongnize that the A-cell contains a function handle instead of numeric matrivuss (and thus not being indicative of the state vector dimension):

%the fragment (lines 700--707) from bnlssm.m that calls the error msg and halts further execution :

if iscellA; At = A{1}; else; At = A; end
numStates0 = size(At,2);
if isscalar(mean0)
    mean0 = mean0(ones(numStates0,1),:);
elseif ~isempty(mean0)
    mean0 = mean0(:);
    validateattributes(mean0,{'numeric'},{'finite','real','numel',numStates0},callerName,'Mean0');
end

The working example i considered is as follows:

rng(1,"twister") % For reproducibility

T = 100;

thetatrue = [0.8; .9; -0.75; -.9; 0.3; 0.2; 0.1];

MdlSim = varm(AR={diag(thetatrue([1 3]))},Covariance=diag(thetatrue(5:6).^2), ...

Constant=thetatrue([2 4]));

XSim = simulate(MdlSim,T);

ysim = log(sum(exp(XSim - mean(XSim)),2)) + thetatrue(7)*randn(T,1);

Mdl = bnlssm(@(params)paramMap(params, otherData), ...

@(params)priorDistribution(params, otherData));

[xhat, SmStats] = smooth(Mdl,ysim,thetatrue);

function logprior = priorDistribution(Params, otherData)

% Prior means for the parameters

prior_means = otherData.prior_means;

% Prior standard deviations for the first 4 parameters

prior_std_devs = otherData.prior_sd;

% Degrees of freedom for the chi-squared distribution (last 3 parameters)

chi2_df =otherData.chi2_df;

% Check parameter constraints

paramconstraints = [(abs(Params([1 3])) < 1) (abs(Params([2 4])) > .2) (abs(Params([1 3])) < 3) (Params(5:7) > 0)];

% keyboard

if all(paramconstraints)

% Compute log-prior for the first 4 parameters (normal distribution)

logprior_normal = sum(log(normpdf(Params(1:4), prior_means(1:4), prior_std_devs)));

% Compute log-prior for the last 3 parameters (chi-squared distribution)

logprior_chi2 = sum(log(chi2pdf(1 ./ Params(5:7), chi2_df)));

% Total log-prior is the sum of log-priors for normal and chi-squared parts

logprior = logprior_normal + logprior_chi2;

else

% If parameter constraints are violated, assign a very low log-prior

logprior = -Inf;

end

end

function [A,B,C,D,Mean0,Cov0,StateType] = paramMap(Params, otherData)

T = otherData.T;

% A = @(x)blkdiag([Params(1) Params(2); 0 1],[Params(3) Params(4); 0 1])*x;

B = [Params(5) 0; 0 0; 0 Params(6); 0 0];

D = Params(7);

C = @(x)log(exp(x(1)-Params(2)/(1-Params(1))) + ...

exp(x(3)-Params(4)/(1-Params(3))));

A = cell(T,1);

for t=1:T

A{t} = @(x)blkdiag([Params(1) Params(2); 0 1],[Params(3) Params(4); 0 1])*x;

end

Mean0 = [Params(2)/(1-Params(1)); 1; Params(4)/(1-Params(3)); 1];

Cov0 = diag([Params(5)^2/(1-Params(1)^2) 0 Params(6)^2/(1-Params(3)^2) 0]);

StateType = [0 1 0 1]; % Stationary state and constant 1 processes

end

the example works fine if 'A' is defined as constant coef. matrix, i.e.:

function [A,B,C,D,Mean0,Cov0,StateType] = paramMap(Params, otherData)
    T = otherData.T;
    A = @(x)blkdiag([Params(1) Params(2); 0 1],[Params(3) Params(4); 0 1])*x; 
    B = [Params(5) 0; 0 0; 0 Params(6); 0 0];
    D = Params(7);
    C = @(x)log(exp(x(1)-Params(2)/(1-Params(1))) + ...
        exp(x(3)-Params(4)/(1-Params(3))));
    Mean0 = [Params(2)/(1-Params(1)); 1; Params(4)/(1-Params(3)); 1];         
    Cov0 = diag([Params(5)^2/(1-Params(1)^2) 0 Params(6)^2/(1-Params(3)^2) 0]);          
    StateType = [0 1 0 1];     % Stationary state and constant 1 processes
end

The same error occurs if i feed as cell 'B' or 'D', which do not feature a function handle for x, e.g.:

function [A,B,C,D,Mean0,Cov0,StateType] = paramMap(Params, otherData)
    T = otherData.T;
    A = @(x)blkdiag([Params(1) Params(2); 0 1],[Params(3) Params(4); 0 1])*x; 
    B = [Params(5) 0; 0 0; 0 Params(6); 0 0];
    D = Params(7);
    C = @(x)log(exp(x(1)-Params(2)/(1-Params(1))) + ...
        exp(x(3)-Params(4)/(1-Params(3))));
    D = cell(T,1);
    for t=1:T
        D{t} = Params(7);
    end
    Mean0 = [Params(2)/(1-Params(1)); 1; Params(4)/(1-Params(3)); 1];         
    Cov0 = diag([Params(5)^2/(1-Params(1)^2) 0 Params(6)^2/(1-Params(3)^2) 0]);          
    StateType = [0 1 0 1];     % Stationary state and constant 1 processes
end

or if feed only only selected elements of At, as function handles.

Please advise me on a work-around that error message and/or give a working example of 'bnlssm' object's application that does feature A_t(x_{t-1}), wheby the coefficients are differnent across periods while also being a non-linear functions of the previos period's state vector elements.