How to use the DeflatedData option in bnlssm package

I hope to use the Bayesian Non-linear state space model in the matlab which names bnlssm. I set a simple model to learn it. I want to add predictors variables, that is, use the DeflatedData option, but it keeps failing. I would like to ask how to use this option. Below is my code.

data = readtable('C:\Users\Xiaoxuan\SSE.csv');
head(data);
SSE = data.close;
dts = data.time;
dts = datetime(dts, 'InputFormat','MM/dd/yyyy');
T = numel(SSE);

retsp500 = price2ret(SSE);
y = retsp500 - mean(retsp500);
retsse = price2ret(SSE);
retdts = dts(2:end);

Z = [data.high];

figure
plot(dts,SSE)
title("China SSE Closing Prices")
ylabel("Closing Price")



PriorMdl = bnlssm(@(theta)paramMap(theta,y,Z),@flatLogPrior,ObservationForm="distribution", ...
    Multipoint=["A" "LogY"]);



theta0 = [0.2 0 0.5 0.7 0 1 1];
lower = [-1; -Inf; 0; 0; -Inf; -Inf; -Inf];
upper = [1; Inf; 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,DeflatedY] = paramMap(theta,y,Z)
    A = @(x) theta(2) + theta(1).*x;
    B = theta(3);
    LogY = @(y,x) -0.5*log(2*pi) - x - 0.5*log(theta(4)) - 0.5*(y-theta(5))*(y-theta(5)).*exp(-2.*x)/theta(4);
    Mean0 = theta(2)./(1 - theta(1));          
    Cov0 = (theta(3)^2)./(1 - theta(1)^2);
    StateType = 0;
    DeflatedY = y - Z*[theta(6); theta(7)];
end



function logprior = flatLogPrior(theta)
% flatLogPrior computes the log of the flat prior density for the three
% variables in theta. Log probabilities for parameters outside the
% parameter space are -Inf.
    paramcon = zeros(numel(theta),1);
    % theta(1) is the lag 1 term in a stationary AR(1) model. The
    % value needs to be within the unit circle.
    paramcon(1) = abs(theta(1)) >= 1 - eps;

    % alpha must be finite
    paramcon(2) = ~isfinite(theta(2)); 

    % Standard deviation of the state disturbance theta(3) must be positive.
    paramcon(3) = theta(3) <= eps;

    % Standard deviation of the state disturbance theta(3) must be positive.
    paramcon(4) = theta(4) <= eps;

    % alpha must be finite
    paramcon(5) = ~isfinite(theta(5)); 

    % alpha must be finite
    paramcon(6) = ~isfinite(theta(6)); 

    % alpha must be finite
    paramcon(7) = ~isfinite(theta(7)); 

    if sum(paramcon) > 0
        logprior = -Inf;
    else
        logprior = 0;   % Prior density is proportional to 1 for all values
                        % in the parameter space.
    end
end

[Ronit]

Hi @Yuanqing, Please share the data file to look into this further.

[Yuanqing]

Hi @Ronit,

Thanks for your reply. The "SSE.csv" in the attachment is a test data with few observations. If you can't open it, please let me know.

I chose the close price as the observed variable y. At the same time, in order to learn how to use DeflatedData option to add explanatory variables, I simply chose the high price variable for analysis, names Z. And the final prior distribution is also set arbitrarily.

Answer 1

[Ronit]

Hi @Yuanqing,

I understand that you are working with a Bayesian Non-linear State Space Model (bnlssm) and want to incorporate a deflation of your observed data using a predictor variable. Please consider the following recommendations to correctly implement the code:

1. Data Alignment: Make sure that the length of variable "Z" matches the length of variable "y".
2. Deflation: In the "paramMap" function, perform deflation in accordance with the dimension of "Z".
3. Parameter Mapping: Make sure that the number of parameters in "theta" corresponds to the operations being performed in "paramMap". You have "theta(6)" and "theta(7)" for deflation, but if "Z" is just one variable, you might only need "theta(6)".

Please review the following code, which accommodates the above mentioned changes with appropriate comments:

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

% Ensure Z is the same length as y
Z = data.high(2:end);

figure
plot(dts, SSE)
title("China SSE Closing Prices")
ylabel("Closing Price")

PriorMdl = bnlssm(@(theta)paramMap(theta, y, 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, DeflatedY] = paramMap(theta, y, Z)
    A = @(x) theta(2) + theta(1).*x;
    B = theta(3);
    LogY = @(y, x) -0.5*log(2*pi) - x - 0.5*log(theta(4)) - 0.5*(y-theta(5))*(y-theta(5)).*exp(-2.*x)/theta(4);
    Mean0 = theta(2)./(1 - theta(1));          
    Cov0 = (theta(3)^2)./(1 - theta(1)^2);
    StateType = 0;

    % Adjusted for single predictor variable
    DeflatedY = y - Z * theta(6);
end

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

Presented below is the output:

         Optimization and Tuning        
      | Params0  Optimized  ProposalStd 
----------------------------------------
 c(1) |  0.2000    0.7038      0.2511   
 c(2) |   0        0.1582      0.7160   
 c(3) |  0.5000    1.1916      0.5439   
 c(4) |  0.7000    0.0000      0.0000   
 c(5) |   0        0.0691      0.0145   
 c(6) |   1       -0.0006      0.0001   
 
      Posterior Distributions of Parameters     
      |   Mean     Std   Quantile05  Quantile95 
------------------------------------------------
 c(1) |  0.5613  0.2570     0.0735      0.8938  
 c(2) | -0.6103  0.5168    -1.6142      0.0219  
 c(3) |  0.8223  0.2907     0.4307      1.3647  
 c(4) |  0.0005  0.0003     0.0000      0.0008  
 c(5) |  0.2299  0.0409     0.1711      0.3021  
 c(6) | -0.0018  0.0003    -0.0023     -0.0014  

     Posterior Distributions of Final States    
      |   Mean     Std   Quantile05  Quantile95 
------------------------------------------------
 x(1) | -1.8477  1.0618    -3.4045     -0.0935  
Proposal acceptance rate = 12.06%

I hope it helps your query!

[Yuanqing]

Hi @Ronit,

Thank you so much for your detailed answer! Your suggestions were spot on, and they have completely resolved my issue. The code works perfectly now.