如何在R中使用hmm TMB提前一步预测观察到的状态

目前，hmmTMB中没有内置的预测功能(截至2024年2月).然而，预测分布可以从拟合的hmm TMB模型计算出来.

预测分布是状态依赖分布的混合，其中权重对应于处于不同状态的概率.因此，您可以使用以下工作流:

1. 获取最后观察到的数据行的状态分布(即，在时间n)；让我们称其为u
2. 得到时间n+h的状态分布为u %*% (tpm %^% h)，其中tpm是状态过程的(一步)转移概率矩阵；也就是说，我们将u乘以h步转移概率矩阵
3. 将预测分布作为估计的状态相关分布的混合，按上一步中找到的状态概率进行加权

例如，在步骤3中，您可以计算预测平均值(作为状态相关平均值的加权和)，或预测分布的概率密度函数，如下面的代码所示.

预测均值与时间

黑点表示训练集最后一次观测后60天的预测分布的平均值.

预测分配

这些线条显示了训练集最后一次观测后60天的预测分布的概率密度函数.正如您所看到的，它很可能一开始处于状态2，因此相应的状态依赖分布具有更高的权重.一段时间后，随着状态分布接近平稳(长期)分布，预测分布停止变化.

代码

# Load packages
library(scico)
library(expm)

# Grid of oil prices to plot forecast distributions
grid <- seq(min(training_set$OILPRICE), 
            max(training_set$OILPRICE), 
            length = 100)

# Get state distribution at last (training) observation
sp <- hmm1$state_probs()
u <- sp[nrow(sp),]

# Get model parameters
tpm <- hid1$tpm()[,,1]
obspar <- obs1$par()[,,1]

# Loop over forecast times (over 60 days here)
h <- 1:60
mix <- list()
mix_mean <- rep(NA, length = length(h))
for(i in seq_along(h)) {
    # State distribution at time n + h
    dist_h <- u %*% (tpm %^% h[i])
    
    # 预测分配 at time n + h
    mix[[i]] <- 
        dist_h[1] * dnorm(x = grid, 
                          mean = obspar["OILPRICE.mean", "state 1"], 
                          sd = obspar["OILPRICE.sd", "state 1"]) +
        dist_h[2] * dnorm(x = grid, 
                          mean = obspar["OILPRICE.mean", "state 2"], 
                          sd = obspar["OILPRICE.sd", "state 2"])
    
    # Mean of forecast distribution
    mix_mean[i] <- sum(dist_h * obspar["OILPRICE.mean",])
}

# Plot the forecast mean against time
df_mean <- data.frame(day = seq.Date(from = as.Date("2020-10-20"), 
                                     by = "day", length = length(h)),
                      mean = mix_mean)

ggplot(data_plot, aes(day, OILPRICE, col = state)) +
    geom_point() +
    geom_point(aes(y = mean, col = NULL), data = df_mean) +
    scale_color_manual(values = pal)

# Plot the forecast distributions for h = 1, 2, ..., 60
df_mix <- data.frame(h = rep(h, each = 100),
                 price = grid,
                 mix = unlist(mix))

ggplot(df_mix, aes(price, mix, col = h, group = h)) +
    geom_line() +
    scale_color_scico(name = "time") +
    labs(y = "forecast distribution")

编辑

请注意，我使用了一个略有不同的模拟数据集，该数据集是从2状态隐马尔可夫模型生成的，以创建我们期望在真实数据中找到的那种模式.数据可以模拟如下:

# Set seed for reproducibility
set.seed(123)

# Generate state process for simulated data
markov <- rep(1, 366)
for(i in 2:366) {
    switch <- sample(0:1, size = 1, prob = c(0.95, 0.05))
    if(switch) markov[i] <- 3 - markov[i-1]
    else markov[i] <- markov[i-1]
}

# Example dataset preparation (extended for 1 year)
df <- data.frame(
    day = seq.Date(from = as.Date("2020-01-01"), 
                   to = as.Date("2020-12-31"), 
                   by = "day"),
    
    # Generate observations conditionally on state process
    OILPRICE = rnorm(366, c(100, 140)[markov], 10)
)