How to use the DeflatedData option in bnlssm package

4 次查看（过去 30 天）
Yuanqing
Yuanqing 2024-9-26
评论： Yuanqing 2024-9-27
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
  2 个评论
Yuanqing
Yuanqing 2024-9-27
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Ronit
Ronit 2024-9-27
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!
  1 个评论
Yuanqing
Yuanqing 2024-9-27
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Weather and Atmospheric Science 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by